Markov  chain-based  policies  for  multi-stage  stochastic  integer  linear  programming

Markov Chain-based Policies for Multi-stage Stochastic Integer Linear Programming with an Application to Disaster Relief Logistics

Margarita P. Castro
Department of Industrial and Systems Engineering, Pontificia Universidad Católica de Chile, Santiago 7820436, Chile<EMAIL_ADDRESS>Merve Bodur
Department of Mechanical and Industrial Engineering, University of Toronto, Toronto, Ontario M5S 3GH, Canada<EMAIL_ADDRESS>Yongjia Song
Department of Industrial Engineering, Clemson University, Clemson, South Carolina 29631<EMAIL_ADDRESS>
We introduce a novel aggregation framework to address multi-stage stochastic programs with mixed-integer state variables and continuous local variables (MSILPs). Our aggregation framework imposes additional structure to the integer state variables by leveraging the information of the underlying stochastic process, which is modeled as a Markov chain (MC). We present a novel branch-and-cut algorithm integrated with stochastic dual dynamic programming as an exact solution method to the aggregated MSILP, which can also be used in an approximation form to obtain dual bounds and implementable feasible solutions. Moreover, we apply two-stage linear decision rule (2SLDR) approximations and propose MC-based variants to obtain high-quality decision policies with significantly reduced computational effort. We test the proposed methodologies in a novel MSILP model for hurricane disaster relief logistics planning. Our empirical evaluation compares the effectiveness of the various proposed approaches and analyzes the trade-offs between policy flexibility, solution quality, and computational effort. Specifically, the 2SLDR approximation yields provable near-optimal solutions for our test instances supported by the proposed bounding procedure.

Multi-stage stochastic programming, Markov chain, linear decision rules, stochastic dual dynamic programming, disaster relief logistics

§ INTRODUCTION

In this paper, we focus on multi-stage stochastic integer linear programming (MSILP) problems, a class of mathematical programming models for sequential decision-making under uncertainty. These problems consider decisions at each stage of the sequential process, which consist of local variables—variables that only participate in a single stage locally—and state variables—variables that link multiple stages together. In particular, we consider MSILP problems with continuous local variables and  mixed-integer state variables, which arise in many applications such as hydro-power scheduling [Hjelmeland et al., 2018], unit commitment [Zou et al., 2018], and disaster relief logistics planning (see Section <ref>).

Despite the vast applicability of MSILP, these problems are computationally prohibitive to solve due to the non-convexity of the feasibility set caused by the integer state variables. Existing solution methods utilize the stochastic dual dynamic programming (SDDP) algorithm with convexification procedures (see, e.g., [Zou et al., 2019] and [Ahmed et al., 2020]), which require strong assumptions about the feasibility set (e.g., binary state variables and Lipschitz continuous cost-to-go functions).

We propose an alternative approach to handle such MSILP problems. The main idea is to construct a partially extended reformulation by transforming the original problem into one that has integer variables only in the first stage. By doing so, we can employ decomposition algorithms, such as the branch-and-cut (B&C) algorithm, to decompose the problem into a master problem that corresponds to a first-stage mixed integer linear program, and a subproblem defined for the remaining stages, which is a multi-stage stochastic linear program that can be handled by stage-wise decomposition algorithms such as the SDDP algorithm. However, this reformulation may lead to exponentially many first-stage integer variables, which is computationally challenging to handle.

To alleviate the main challenge of the partially extended reformulation, we present a novel aggregation framework that imposes structures to the reformulation by aggregating the integer state variables while potentially preserving high-quality solutions. This framework leverages the structure of the underlying stochastic process to impose a restriction on the problem at hand, leading to an informative decision policy. In particular, we consider the case where the underlying stochastic process is modeled as a Markov chain (MC) and present several transformations based on the MC state information that lead to a wide range of aggregated MSILP approximations.

In order to solve the resulting aggregated MSILP models, we propose a novel B&C framework integrated with a variant of the SDDP algorithm. The framework decomposes the problem into a master problem and a set of subproblems, where the master problem is concerned with the first-stage decisions including all the integer state variables, and the subproblems deal with the remaining stages, which involve continuous state and local variables. This exact solution method can also be used to compute relaxation bounds with relaxed termination criteria.

While the proposed B&C framework guarantees an optimal solution for the aggregated MSILP, solving a multi-stage stochastic program on every node of the branch-and-bound tree can be computationally challenging. Therefore, we consider a more tractable alternative that approximates the aggregated MSILP as a two-stage stochastic program using two-stage linear decision rules (2SLDRs) [Bodur  Luedtke, 2022]. We present 2SLDR alternatives from the literature and propose a novel variant that leverages the MC structure of the underlying stochastic process. The resulting problem can be solved via Benders decomposition or any of its enhanced versions.

We perform a case study for the proposed aggregation framework and solution methods on a novel hurricane disaster relief logistic planning problem (HDR), where the stochastic process is given by an MC that models the evolution of the hurricane. We show how to create several aggregated MSILPs and 2SLDR approximations for this application. Our numerical experiments compare these alternatives, find suitable transformations and 2SLDR variants to obtain high-quality solutions with provable bounds for this application, and generate insights about the trade-off between solution quality and computational effort.

The remainder of the paper is organized as follows. Section <ref> presents our aggregation framework and a literature review on MSILP methodologies. Section <ref> introduces our novel B&C framework integrated with the SDDP algorithm and describes the 2SLDR scheme with our proposed MC-based variant. Section <ref> presents  a case study of the proposed methodologies for HDR. Section <ref> presents the numerical experiment results and we end with some concluding remarks in Section <ref>.

* This work focuses on multi-stage stochastic integer linear programming (MSILP) problems with mixed-integer state variables and continuous local/recourse variables.

* These types of problems are very challenging, arise in many applications, and there are limited procedures to solve them. One of the main challenges with these types of problems is the presence of integer variables since the feasibility set is nonconvex.

* If  we only have integer variables in the first stage, we can handle this problem. We can decompose the problem into two phases. The first phase is a master problem and the second phase is a MSP with only continuous variables. It can also be solved by means of approximation using 2S-LDR.

* There are few techniques to tackle these problems when we have integer variables in other stages. However, they are limited to special cases (e.g., binary variables, Lipschitz continuous). These techniques are variants of SDDP to handle integer variables.

* With these in mind, we would like to transform our problem into the case when we only have integer variables in the first stage. This can be achieved by creating scenario copies of the integer variables and adding them to the first-stage problem along with nonanticipativity constraints. If the uncertainty model is in the form of a scenario tree, an alternative would be to create node-based copies of those variables instead, that are nonanticipative by definition. We call these variables first-stage reformulation variables.
Such a reformulation usually requires an exponential number of new variables (even when we have stage-wise independence) since they depend on the history of the stochastic process. As such, it would be very challenging to solve.

* One way to alleviate this difficulty is to remove the dependence of the  first-stage reformulation variables to the history, by restricting them in a certain way to reduce their number. For instance, we can make them all the same regardless of their history if they belong to the same stage (i.e., stage-based restriction). Note that this is different to the standard two-stage restriction since we leave the continuous stage variables the same (reference from Shabbir). (find some references for the integer transformation—very likely that there are not).

* The stage-based restriction would usually lead to poor quality solutions.
Thus, we extend this idea and propose MC-based restrictions for  problems where the uncertainty is modeled as a Markov chain (see Section X.X for generalizations and stage-wise independent case).
That is, we try to leverage the structure of the underlying stochastic process to create high-quality policies.
We restrict the first-stage reformulation variables to be a function of the “Markov state" of their associated stage (or their limited history), to get a good balance of the difficulty of solving the reformulation and the quality of the policy.
Thus, the resulting restriction will capture the structure of the stochastic process and allow us to have a manageable number of integer first-stage variables and only continuous variables in the other stages.

* In order to solve the obtained restriction, where we have integer variables only in the first stage, we propose a branch-and-cut (B&C) algorithm. The master problem only considers the first-stage variables. The subproblem is an MSP with only continuous variables that can be solved, for instance, with an SDDP algorithm. To the best of our knowledge, B&C algorithms where the node relaxation is solved by SDDP
haven't been applied before to MSILP with only integer first-stage variables.

* Solving an MSP on every node of the branch-and-bound tree can be very challenging. Therefore, we also consider a more tractable alternative that transforms our restriction problem into a two-stage SP via LDR, by restricting the continuous state variables to be a linear function of the history of the stochastic process. We consider some of the two-stage LDR presented in the literature and propose new ones using similar ideas to create policies for the integer state variables. The resulting problem can be solved via Benders decomposition and any of its enhancements proposed in the literature.

* We test our procedure in a disaster relief problem for hurricanes, where the stochastic process is defined by a Markov chain. We show how to create MC-based policies and LDRs for this application. We compare them again stage-based policies and LDR, and the optimal policy solved by an extensive formulation.

* Then, we have an experimental section focusing on the hurricane application.

* Finish with the conclusions.

§ AGGREGATION FRAMEWORK FOR INTEGER STATE VARIABLES

This section describes the proposed aggregation framework to generate informative policies for a class of MSILP problems. As previously mentioned, the main idea is to aggregate integer state variables in the MSILP by leveraging the structure of the underlying stochastic process.

§.§ Problem Formulation for MSILP

We consider a class of MSILP models with continuous local variables and mixed-integer state variables. Although integer variables may also appear as local variables in a generic MSILP, the restriction here (local variables being continuous only) does not limit the model from being applicable to a broad range of problems that arise in real-world applications. For example, a typical multi-period stochastic optimization problem with local production and distribution constraints can be represented with this model, such as the one we study in Section <ref>, motivated by the application of disaster relief logistics planning. Other possible applications include hydro-power scheduling [Hjelmeland et al., 2018] and unit commitment [Zou et al., 2018] problems.

We assume that the underlying stochastic process has finite support and, as such, it can be represented by a scenario tree $\tree$. This assumption is not restrictive since we can construct an approximate scenario tree if the stochastic process has continuous support by means of sampling, a common practice in the literature [Shapiro et al., 2021]. Let $\periods$ be the number of stages and let $\nodes$ be the set of nodes associated with the scenario tree $\tree$. The set of nodes in each stage $t \in \horizon:=\{1,...,\periods\}$ is given by $\staget$. The root node is denoted by $\rootnode$ and it is the only node in the first stage, i.e., $\stage_1=\{\rootnode\}$. Each node $n \in \staget$ with $t>1$ has a unique parent in $\stage_{t-1}$, denoted by $\parent(n)$. Thus, there is a unique path from $\rootnode$ to any node $n\in \staget$ in stage $t>1$, and we let $\patht(n)$ represent the set of nodes on this path (including $\rootnode$ and $n$). For each non-leaf node $n$ (i.e., $n\in \staget$ for $t \in [\periods-1]$), we let $\childs(n)$ be the set of children of $n$, that is, the nodes whose parent is $n$. We denote the probability that node $n$ occurs by $\prob_n$, and it is clear that we have $\prob_n>0$ for all $n\in \nodes$ and $\sum_{n \in \staget} \prob_n = 1$ for all $t\in \horizon$. The probability of transitioning from node $n \in \staget$ for $t\in [\periods-1]$ to a child node $n' \in \childs(n)$ can then be written as $\pprob_{nn'}:= \prob_{n'}/\prob_n$.

Next, we present a generic problem formulation for the class of MSILP problems that we consider in (<ref>), which is often referred to as the nested formulation. For each node $n\in \nodes$, let $\xn \in \R^k$ and $\zn \in \Z^\ell$ be the set of continuous and integer state variables, respectively, and $\yn\in \R^r$ be the continuous local variables. The objective coefficient vectors associated with $\xn, \yn$, and $\zn$ are denoted by $d_n, h_n$, and $c_n$, respectively. The overall objective is to minimize the sum of the local cost and the expected future cost at the root node:

\begin{align}
Q_{\rootnode} = \min\; & c_\rootnode^\top \zroot + d_\rootnode^\top \xroot + h_\rootnode^\top \yroot + \sum_{n \in \childs(\rootnode)}\pprob_{\rootnode n} Q_n(\xroot,\zroot ) , \tag{$P$} \label{mod:scenarioTree} \\
\text{s.t.}\ &
H_{\rootnode} \zroot \geq g_{\rootnode}, \;
J_{\rootnode} \xroot \geq f_{\rootnode}, \;
C_{\rootnode} \xroot + D_{\rootnode} \zroot + E_{\rootnode} \yroot  \geq  b_{\rootnode}, \label{eq:st_root}\\
&\xroot \in \R^k,\; \yroot \in \R^r,\; \zroot \in \Z^\ell, \nonumber
\end{align}

where the cost-to-go function $Q_n(\cdot, \cdot)$ associated with each non-root node $n\in \nodes\setminus \{\rootnode\}$ is defined in the following nested fashion:

\begin{align}
Q_n(\xnparent, \znparent) = \min\; &  c_n^\top \zn + d_n^\top \xn + h_n^\top \yn + \sum_{n' \in \childs(n)}\pprob_{nn'} Q_{n'}(\xn,\zn),\nonumber \\
\text{s.t.}\ & H_n \zn \geq  G_n \znparent + g_n, \label{eq:st_onlyz} \\
& J_n \xn \geq  F_n \xnparent + f_n, \label{eq:st_onlyx} \\
& C_n \xn + D_n \zn + E_n \yn  \geq A_n \xnparent + B_n  \znparent + b_n, \label{eq:st_all} \\
&\xn \in \R^k,\; \yn \in \R^r,\; \zn \in \Z^\ell. \nonumber
\end{align}

The cost-to-go function $Q_n(\cdot, \cdot)$ associated with a leaf node $n\in \stage_\periods$ only involves the cost incurred in the terminal stage $\periods$ since its set of children nodes is empty, i.e., $\childs(n) = \emptyset$. Our proposed formulation considers, separately, one set of constraints (<ref>) for integer state variables, one set of constraints (<ref>) for continuous state variables, and one set of constraints (<ref>) that link state and local variables. Although alternative formulations are possible, this form accommodates the presentation of the proposed aggregation framework, which will be described in Section <ref>.

Problem (<ref>) is notoriously challenging to solve for two reasons. The first challenge is the nature of multi-stage stochastic programs where the underlying stochastic process is modeled as a scenario tree of an exponentially large size. Thus, it is impractical to solve these large-scale problems directly and decomposition methods can be deemed necessary. The second challenge is the existence of integer decision variables, which makes the expected cost-to-go functions (i.e., the value functions that represent the expected future cost defined at each node of the scenario tree) to be nonconvex in general. The following section reviews existing methods in the literature to handle (<ref>) and similar problems, focusing on how existing approaches in the literature address these challenges.

§.§ Literature Review

In this section, we discuss recent advances in the literature that address MSILPs like formulation (<ref>). For MSILPs defined with a large-scale scenario tree, decomposition algorithms are necessary to obtain high-quality solutions in a computationally viable fashion.

Nested Benders decomposition [Birge, 1985] and its sampling variant under the assumption of stage-wise independence, SDDP [Pereira  Pinto, 1991], are classical decomposition algorithms for multi-stage stochastic linear programs.
Recently, [Zou et al., 2019] develop the SDDiP algorithm for MSILPs with pure binary state variables, where the nonconvex expected cost-to-go functions are lower approximated by a piecewise linear convex envelope constructed by the so-called Lagrangian cuts. These cuts are guaranteed to be exact at points where the associated state variables are binary, leading to the exactness and finite convergence of the  algorithm. The problem becomes more challenging when the set of state variables includes general mixed-integer decision variables, which is the case for formulation (<ref>). One viable approach is to perform the so-called binarization procedure to reformulate general integer state variables and approximate continuous state variables via binary variables, respectively, and then apply the SDDiP algorithm [Zou et al., 2018]. This approach may suffer from a large number of binary state variables as a result of the binarization.

Alternatively, when the expected cost-to-go functions are assumed to be Lipschitz continuous, [Ahmed et al., 2020] propose nonlinear cuts to exactly approximate the nonconvex expected cost-to-go functions via the augmented Lagrangian method. More recently, [Zhang  Sun, 2019] propose and study a unified theoretical framework for general multistage nonconvex stochastic mixed integer nonlinear programming problems where the expected cost-to-go functions are not necessarily Lipschitz continuous. Specifically, they propose regularized expected cost-to-go functions that can be exactly approximated by the generalized conjugacy cuts, generalizing the results of [Ahmed et al., 2020]. However, despite these theoretical convergence properties, the nonlinear cut generation from augmented Lagrangian and the implementation of a variant of the nested Benders decomposition algorithm that incorporates these cuts can be computationally challenging.

More recently, [Füllner  Rebennack, 2022] extend the ideas of binarization procedure in [Zou et al., 2019] and the regularization procedure in [Zhang  Sun, 2019] to develop the so-called nonconvex nested Benders decomposition algorithm for generic multi-stage mixed integer nonconvex nonlinear program. Despite its generality and successful implementation in deterministic multi-stage MINLP problems, it is unclear, to the best of our knowledge, whether or not applying this approach to MSILPs with a large-scale scenario tree is computationally feasible.

In sum, all the works mentioned above aim to directly address the nonconvex (and sometimes non-Lipschitz) expected cost-to-go functions defined at each node of the scenario tree by developing exact lower-bounding techniques. In contrast, we circumvent the challenge of approximating the nonconvex expected cost-to-go functions by relocating all integer state variables (defined in the node subproblems) to the first stage (defined in the root node), so that the resulting cost-to-go functions defined at non-root nodes are convex and can be approximated (exactly) by a decomposition scheme (e.g., nested Benders or SDDP). We next discuss our proposed reformulation in detail.

§.§ A Partially Extended Formulation for MSILP

We propose a reformulation for MSILP in the form of a partially extended formulation, in which all integer decisions can be considered as the first-stage variables. The advantage of this reformulation is two-fold. First, the resulting cost-to-go functions at each node $n \in \nodesnoroot$ have only continuous variables (including local and continuous state variables), making it amenable to apply standard Benders-type cutting-plane approximation to these convex functions. Second, since integer variables only appear in the first-stage problem, they can be handled with a B&C procedure (see Section <ref>). Specifically, the proposed partially extended (nested) formulation is given by:

\begin{align*}
\Qref_{\rootnode} = \min\; &
\sum_{n\in\nodes}\prob_nc_n^\top \zn
+ d_\rootnode^\top \xroot + h_\rootnode^\top \yroot + \sum_{n \in \childs(\rootnode)} \pprob_{\rootnode n} \Qref_n(\xroot,z), \tag{$P^R$} \label{mod:reformulation} \\
\text{s.t.}\; &
H_{\rootnode} \zroot \geq g_{\rootnode}, \;
J_{\rootnode} \xroot \geq f_{\rootnode}, \;
C_{\rootnode} \xroot + D_{\rootnode} \zroot + E_{\rootnode} \yroot  \geq  b_{\rootnode},  \\
& H_n \zn \geq  G_n \znparent + g_n, & \forall n \in \nodesnoroot,
\end{align*}

where the cost-to-go function $Q^R_n(\cdot, \cdot)$ associated with
$n\in \nodes\setminus \{\rootnode\}$ is defined in a  nested fashion as:

\begin{equation*}
\Qref_n(\xnparent, z) = \min_{\xn\in \R^k, \yn \in \R^r} \left\{d_n^\top \xn + h_n^\top \yn  + \sum_{n' \in \childs(n)} \pprob_{n n'} \Qref_{n'}(\xn, z) \mid \eqref{eq:st_onlyx}-\eqref{eq:st_all}\right\}.
\end{equation*}

The key difference between this reformulation (<ref>) and the original formulation (<ref>) is that all the integer state variables $z = \{\zn\}_{n \in \nodes}$, along with their associated constraint sets (<ref>), are now “moved” to the problem at the root node $\rootnode$, i.e., they are treated as first-stage variables and constraints. As such, the cost-to-go function associated with each non-root node only includes the state and local continuous variables and their associated constraints. Since the continuous state variables remain to be defined in the individual node-based subproblems through the nested form, we call this reformulation a partially extended formulation. The validity of this reformulation is clear because we only relocate the (node-based) integer state variables in order to achieve the desired property of piece-wise linear convex expected cost-to-go functions, which accommodates the application of stage-wise decomposition methods such as SDDP.

§.§ An Aggregation Framework via MC-based Structural Policies

While the proposed partially extended formulation (<ref>) has piece-wise linear convex expected cost-to-go functions, it may lead to exponentially many integer variables in the first stage, which is computationally challenging to handle. We propose a novel approximation framework that imposes structures to the formulation by aggregating the integer state variables to alleviate this challenge while maintaining high-quality solutions. The idea is to leverage the structure of the underlying stochastic process to impose a restriction to the problem and obtain an informative decision policy.

In this paper, we consider the case where the underlying stochastic process is modeled as a Markov chain (MC). We choose this specific structure for two reasons. First, many stochastic processes in real-world applications can be modeled or approximated using MCs [Löhndorf  Shapiro, 2019], including disaster relief logistics planning problems that we focus on in our numerical experiments (see Section <ref>). Second, it facilitates the exposition of our aggregation framework and allows us to construct a wide range of structured policies. We discuss alternative structural policies beyond those based on MC models in Remark <ref>.

To start, we let $\markov$ be the set of MC states. Each MC state $m\in \markov$ corresponds to an $\dimmarkov$-dimensional vector, i.e., $m \in \R^\dimmarkov$. Let $\markov_t$ represent the set of MC states that are reachable in stage $t\in \horizon$, starting from a given initial state. We use notation $m^t=(m_1,...,m_t)$, where $m_{t'} \in \markov_{t'}$ for $t'\in [t]$, to represent a sequence of admissible MC states, that is, the probability to reach $m_{t'+1}$ given $m_{t'}$ is positive for all $t' \in [t-1]$.

We now describe the connection between the MC states and the scenario tree model used in the MSILP formulation (<ref>). Each node $n\in \staget$ ($t\geq 1$) of the scenario tree corresponds to an MC state, $m_t(n) \in \markov_t$, and is uniquely determined by the trajectory of the stochastic process from the root node (initial MC state) to MC state $m_t(n)$, that is, the sequence of MC states $m^t(n) = (m_1(\rootnode), ..., m_{t-1}(a(n)),m_{t}(n))$. To simplify the exposition, we write $m^t$ instead of $m^t(n)$ and $m_t$ instead of $m_t(n)$ when we refer to any node $n \in \staget$.

As stated above, the MC representation of the stochastic process allows us to parameterize each node in the scenario tree by means of their corresponding stage and sequence of MC states, i.e., for each $n \in \staget$, we have $n \equiv (t, m^t)$. We can then rewrite the integer state variables $\zn$ associated with node $n$ as $\zn \equiv \ztM$. Based on this parameterization, our aggregation framework relies on a linear transformation that will compress the full MC history up to stage $t$, namely the vector $m^t$ of size $\dimmarkov\cdot t$, into a vector $\maptmarkov m^t$ of size $q_t \leq \dimmarkov\cdot t$, by means of a transformation matrix $\maptmarkov \in  \Z^{ q_t \times \dimmarkov\cdot t}$. That is, it parametrizes the integer state variables of stage $t$ by the limited information $\maptmarkov m^t$:

\begin{equation} \label{eq:generalRestricition}
\zn \equiv \ztM \rightarrow \ztphiM  \qquad  \qquad \forall n \in \staget, t \in \horizon.
\end{equation}

It is clear that this aggregation significantly reduces the number of integer state variables if $q_t \ll \dimmarkov\cdot t$.

Different transformations can impose different structural properties on the discrete decisions and, thus, lead to policies with different levels of aggregation. Table <ref> summarizes a set of transformations, the resulting aggregated variables, and the conditions under which the original integer state variables associated with nodes $n$ and $n'$ in $\stage_t$ correspond to the same (aggregated) variables, i.e., $\zn$ and $\znprime$ have the same value given by $\ztphiM$. We explain each of these policies below.

Name     Transformation matrix $\maptmarkov$     Aggr. var.      Condition for $\zn \equiv \znprime, \ \forall n,n' \in \staget$

$\bzero^{1 \times \dimmarkov \cdot t}$       $\za_{t,0}$      $t(n) = t(n')$

\bzero^{\dimmarkov \times \dimmarkov \cdot (t-1)} \ | \ \identity^{\dimmarkov \times \dimmarkov}\right]$      $\za_{t,m_t} $      $m_t(n) = m_t(n')$

\bzero^{2 \dimmarkov \times \dimmarkov \cdot (t-2)} \ | \ \identity^{2 \dimmarkov \times 2 \dimmarkov}
\right]$       $\za_{t,(m_{t-1},m_t)}$      $ m_t(n) = m_t(n') \ \& \ m_{t-1}(a(n)) = m_{t-1}(a(n'))$

\bzero^{(\bar{\dimmarkov} + \dimmarkov) \times (\dimmarkov \cdot t - \bar{\dimmarkov} - \dimmarkov)} \ | \ \identity^{(\bar{\dimmarkov} + \dimmarkov) \times (\bar{\dimmarkov} + \dimmarkov)}
\right] $     $\za_{t,(\bar{m}_{t-1},m_t)}$     $ m_t(n) = m_t(n') \ \& \ \bar{m}_{t-1}(a(n)) = \bar{m}_{t-1}(a(n')) $

$\identity^{\dimmarkov \cdot t \times  \dimmarkov \cdot t} $     $\za_{t,m^t}$     $ m^t(n) = m^t(n')$

MC-based transformations examples, the resulting aggregated variables (Aggr. var.), and the equivalence of the original variables. $\bzero$ and $\identity$ denote the zero and identity matrices, respectively.

The simplest transformation we consider is what we call a Here-and-Now () transformation, which aggregates
integer state variables associated with nodes that are in the same stage, i.e., they are assumed to take the same values: $\zn \equiv \ztM \rightarrow \za_{t,0}$ for all $n \in \staget$. In other words, the transformation matrix $\maptmarkov = \bzero^{1 \times \dimmarkov \cdot t}$ eliminates all the Markov history information $m^t$, and makes the integer state variables only stage-dependent. The corresponding decision policy as a result of the  transformation will be independent of the realization of the stochastic process, often referred to as a static policy in the literature (see, e.g., [Basciftci et al., 2019]). The  transformation imposes a heavy restriction on the integer state variables and may lead to policies of relatively low quality.

An alternative is to leverage the structure of the MC to create a transformation that leads to more informative policies. For example, the Markovian () transformation aggregates integer state variables associated with nodes that share the same stage and the same MC state: $\zn \equiv \ztM \rightarrow \za_{t,m_t}$ for all $n \in \staget$. That is, the transformation matrix $\maptmarkov$ eliminates all the MC history except for the Markovian state in the current stage $t$, i.e., $\maptmarkov m^t = m_t$.
The intuition behind this transformation is that similar decisions should be made at nodes with the same MC state at the same stage. When it becomes concerning to ignore the history of the MC states in prior stages, one may attempt to retain more information by applying a less restrictive transformation, at a price of a larger number of aggregated variables. One such example is what we refer to as the Double Markovian () transformation, which goes one step beyond  and aggregates integer state variables associated with nodes that share the same MC states in both the current stage and the previous stage. Following the pattern of retaining more MC state information, the extreme case is to consider the entire history of the MC states starting from the initial MC state, referred to as the Full History () transformation. It is clear that the  transformation is equivalent to the original partially extended formulation for MSILP (<ref>), since no integer state variables are effectively aggregated.

Although all the aforementioned MC-based transformations are constructed based on the complete MC state information, we may also consider partial information from MC states. This is advantageous especially if some attributes in the MC state are deemed to be more important than others. Consider a subset of attributes of size $\bar{\dimmarkov} < \dimmarkov$ associated with MC state $m_t$, and assume without loss of generality that these attributes of interest correspond to the last $\bar{\dimmarkov}$ attributes in the MC state vector, i.e., we can represent an MC state as $m_t = (\cdot, \bar{m}_t)\in \R^{\dimmarkov -\bar{\dimmarkov}}\times \R^{\bar{\dimmarkov}}$. Then, we can create a transformation based only on this partial vector $\bar{m}_t$ associated with the MC state $m_t$. As an example, we introduce the  Partial Markovian () transformation that considers  complete MC state information in the current stage but only partial information from the previous stage.

An illustrative example showing nodes with the same aggregated integer state variables for transformations , , and . Nodes with the same color and shape share the same set of integer state variables.

Consider an MC $\markov$ with two possible states (i.e., light 1 and dark 1) and positive transition probabilities for all transitions between states. The first drawing in Figure <ref> depicts a scenario tree built for $\markov$ starting at a light state and representing all possible MC states for $\periods=4$ stages, where the color of each node (light or dark) represents the corresponding MC state. The remaining drawings in Figure <ref> illustrate the effect of aggregation on the nodes of the scenario tree for transformations , , and  (we omit  because in this example the MC state vector has only one attribute). The  tree shows that all nodes in the same stage will share the same set of (aggregated) integer state variables $\za$ (i.e., they have the same color and shape). In contrast,  aggregates integer state variables in nodes with the same MC state and stage. Thus, the  tree shows nodes of two different colors in each stage, which correspond to two separate sets of aggregated integer state variables for each stage (e.g., the transformation for node $n\equiv (t,m^t) = (3,\mcLight{1}\mcDark{1}\mcLight{1})\in \nodes_3$ is $\Phi_3 m^3(\mcLight{1}\mcDark{1}\mcLight{1}) = m_3(\mcLight{1}) = \pentaLightDark$, thus,  $\za_{3,\pentaLightDark}$ represents its associated aggregated variables).
Finally, nodes in each stage of the  tree have up to four different colors, as a result of an aggregation based on the MC states in both the current stage and the previous stage
(e.g., $\za_{3,\pentaDarkDark}$ is associated to node (3,111)  because $\Phi_3 m^3(\mcLight{1}\mcDark{1}\mcDark{1}) = (m_2(\mcDark{1}), m_3(\mcDark{1})) = \pentaDarkDark$).
Overall,  contains four sets of $\za$ variables (i.e., one for each stage),   seven, and   eleven, while the  variant contains 15 sets of $\za$ variables (i.e., one for each node of the scenario tree).

These transformations are just a few examples of a wide range of aggregations that can be constructed within our framework. One key consideration in constructing these aggregations is that there is a clear trade-off between the quality of the underlying decision policy and the corresponding computational effort. Aggregations that include (preserve) more information would naturally lead to better decision policies but also make the resulting problem harder to solve due to the increased number of integer state variables in the first stage. Choosing the right transformation depends on the specific problem type and the available computational resources. We leave the problem of adaptively finding the best transformation as a future research direction.

The main advantage of assuming that the stochastic process is defined by an MC is that most stage-dependent processes can be approximated by an MC using a discretization technique with high accuracy [Bally  Pagès, 2003, Pagès et al., 2004]. Also, the work by [Löhndorf  Shapiro, 2019] shows that such an approximation can lead to high-quality solutions compared to the ones obtained by approximating the stochastic process with autoregressive models in certain applications. We also note that our approach is applicable to stage-wise independent stochastic processes. In these cases, one can approximate the probability distribution of a random variable by partitioning its support into a set of clusters and creating aggregations based on these clusters. For example, considering an MSILP with stochastic demand parameters, we may define three levels of demand (low, medium, and high) and create aggregation schemes based on each demand level.

§.§ Aggregated MSILP

We now present the aggregated MSILP formulation after applying the transformation to the integer state variables $z$ via matrices $\maptmarkov, t\in \horizon$. The new set of integer state variables is given by $\za_t\in \Z^{\ell \cdot q_t}$ for each $t\in \horizon$. We employ a mapping function $\phi_t:\nodes_t\rightarrow \{t\}\times[q_t]$ to relate a node $n\in \staget$, which was originally used to index variables $z$, to the index corresponding to its aggregated variables, i.e., $\phi_t(n) = (t, \maptmarkov m^t)$ for each node $n \in \staget$. Thus, the aggregated MSILP can be written as follows:

\begin{align}
\Qaggre_{\rootnode} = \min\; &  \sum_{n\in \nodes}\prob_n c_n^\top \zphi  + d_\rootnode^\top \xroot + h_\rootnode^\top \yroot + \sum_{n \in \childs(\rootnode)}\pprob_{\rootnode n} \Qaggre_n(\xroot,\za) \tag{$P^A$} \label{mod:aggregated} \\
\text{s.t. }\; &
H_{\rootnode} \zroot \geq g_{\rootnode}, \;
J_{\rootnode} \xroot \geq f_{\rootnode}, \;
C_{\rootnode} \xroot + D_{\rootnode} \zroot + E_{\rootnode} \yroot  \geq  b_{\rootnode},  \label{aggre:root}\\
& H_n \zphi \geq  G_n \zphiparent + g_n, & \forall n \in \nodesnoroot,   \label{aggre:zphi}\\
& \xroot \in \R^k, \yroot \in \R^r, \za \in  \Z^{\ell \cdot \sum_{t \in \horizon} q_t},  \label{aggre:vars}
\end{align}

where the cost-to-go function $\Qaggre_n(\cdot,\cdot)$ associated with
$n\in \nodes\setminus \{\rootnode\}$ is defined in a  nested fashion as:

\begin{align*}
\Qaggre_n(\xnparent, \za) = \min\; & d_n^\top \xn + h_n^\top \yn + \sum_{n' \in \childs(n)}\pprob_{n n'} \Qaggre_{n'}(\xn, \za)\\
\text{s.t.}\; & J_n \xn \geq  F_n \xnparent + f_n, \\
& C_n \xn + D_n \zphi + E_n \yn  \geq A_n \xnparent + B_n  \zphiparent + b_n, \\
& \xn \in \R^k, \yn \in \R^r.
\end{align*}

We note that any feasible solution of (<ref>) is also  feasible to (<ref>) with the same objective value, because the only difference between (<ref>) and (<ref>) is that (<ref>) enforces some integer variables of (<ref>) to take the same values. However, the reverse is not true in general, and hence an optimal solution of (<ref>) might be sub-optimal for (<ref>). This leads to the question of how to find a nontrivial aggregation that results in an optimal policy, which we leave as a future research direction.

Although the aggregated problem (<ref>) potentially has much fewer integer variables than (<ref>), depending on the transformation $\maptmarkov$, it may remain  computationally challenging to solve due to the existence of (aggregated) first-stage integer variables and the continuous state variables defined at each node of the scenario tree. We next introduce solution methodologies to address this issue.

§ SOLUTION METHODOLOGY

We now present two solution methods for the aggregated MSILP formulation (<ref>). Both are decomposition methods implemented as B&C algorithms, taking advantage of the fact that all the integer variables in (<ref>) are in the first stage. The first one is an exact method that employs SDDP as a subroutine. The second method employs 2SLDR to impose additional solution structures that lead to a two-stage approximation to the multi-stage  program. This approximation allows us to employ standard Benders decomposition algorithms to solve the problem, reducing the computational effort without a significant compromise on the solution quality.

§.§ SDDP Integrated B&C Algorithm

We now present an exact algorithm for solving formulation (<ref>) by integrating the SDDP algorithm within a B&C framework. The high-level idea is to decompose the problem into a master problem and a set of subproblems, where the master problem is concerned about the first-stage decisions including all the integer variables of (<ref>), and the subproblems deal with the remaining stages, which involve continuous state and local variables.

The algorithm performs a branch-and-bound (B&B) procedure for the master problem and enters a cut generation procedure every time an integer relaxation solution is encountered (e.g., done through the so-called callback functions in modern IP solvers such as CPLEX). To generate a cut, following the basic idea of Benders decomposition, we need to solve a subproblem associated with each child node of $\rootnode$ in the scenario tree, which deals with the remaining stages of (<ref>) via the SDDP algorithm for a given candidate solution $(\xrootsol, \zasol, \thetasol)$ from the master problem. In particular, we consider the so-called multi-cut version, where a variable $\theta_n$ is introduced for each $n \in \childs(\rootnode)$ to represent the value of the
outer approximation of the cost-to-go function at node $n$. Thus, the master problem is given by:

\begin{align}
\Qsub_{\rootnode} = \min& \sum_{n\in \nodes}\prob_n c_n^\top \zphi  + d_\rootnode^\top \xroot + h_\rootnode^\top \yroot + \sum_{n \in \childs(\rootnode)}\pprob_{\rootnode n} \theta_{n} \\
\text{s.t.}\ & \eqref{aggre:root}-\eqref{aggre:vars}, \nonumber \\
& \theta_n \geq \alpha_n^\top \xroot + \beta_n^\top \za + \gamma_n, \quad \forall  (\alpha_n,\beta_n,\gamma_n)\in \coeffs_n, \ n \in \childs(\rootnode),
\end{align}

where $\coeffs_n$ stores the coefficients associated with all the Benders cuts constructed to approximate $\Qaggre_n(\xroot,\za)$ during the solution procedure. Algorithm <ref> describes the B&C procedure. The procedure receives as input a few parameters used by an SDDP sub-routine to be explained in detail later.

1.2A B&C framework integrated with the SDDP algorithm
[1]
$\algExact, K, \epsilon$
Initialize B&B search for $\Qsub_\rootnode$ (master problem) and the  environment.
Initialize $\coeffs_n :=\emptyset$ for all $n \in \childs(\rootnode)$ and $\algCutAdd := \algFalse$.
Choose a B&B node from its list of open nodes. Solve the associated node relaxation problem.
Current B&B node is infeasible or the relaxation bound is worse than the bound given by an incumbent solution
Prune this B&B node.
Current solution $(\xrootsol, \zasol, \thetasol)$ is integer
%% Cutting plane subroutine
$\algCutAdd := \algFalse$.
$n' \in \childs(\rootnode)$  %% Find cuts parameters for each child subproblem
Get cut coefficients $(\alpha_{n'},\beta_{n'},\gamma_{n'})$ := ($(\xrootsol, \zasol, \thetasol), n',\epsilon,K,\algExact$).
Cut with coefficients $(\alpha_{n'},\beta_{n'},\gamma_{n'})$ is violated by  $(\xrootsol, \zasol, \thetasol)$
$\coeffs_{n'}:= \coeffs_{n'}\cup(\alpha_{n'},\beta_{n'},\gamma_{n'})$, $\algCutAdd:= \algTrue$.
$\algCutAdd = \algFalse$.
Branch on the current  B&B node and update the list of open nodes.
The list of open B&B nodes is empty.

There are a few important points to mention about this procedure. First, we call the SDDP sub-routine for each child of the root node $\rootnode$, in which we perform the SDDP forward pass only considering sample paths that contain that specific child. Also, we look for cuts for all children nodes of the root node $\rootnode$ before re-solving the current B&B node relaxation. An alternative would be to re-solve the B&B node relaxation whenever we find a violated cut, but this option showed worse performance in our preliminary experiments. In addition, we may need to call the SDDP sub-routine multiple times for each candidate integer feasible (relaxation) solution to guarantee that such solution is correctly evaluated (i.e., the $\theta_n$ values for $n \in \childs(\rootnode)$ are accurate representations of their corresponding cost-to-go function values). Lastly, all cuts added to enhance the outer approximation of the cost-to-go functions associated with the SDDP subproblems in the previous iterations are automatically carried over to the current iteration.

We next explain the implementation details of this B&C framework integrated with SDDP. We first describe the SDDP subproblem structure and then present the entire SDDP sub-routine.

§.§.§ SDDP Subproblems.

A key feature of the SDDP algorithm is its capability to leverage the structure of the underlying stochastic process in defining the cost-to-go functions and the corresponding SDDP subproblems. For example, we only need to create one SDDP subproblem per stage in the case of stage-wise independence [Pereira  Pinto, 1991], and one subproblem per stage and MC state in the case when the stage-wise dependence structure is modeled by an MC [Bonnans et al., 2012, Philpott  De Matos, 2012]. In what follows, we describe the SDDP subproblems for formulation (<ref>) based on the MC structure and the employed transformation.

First, we observe that if the first-stage integer variables $\za$ do not exist, then thanks to the Markovian nature of the stochastic process, instead of defining one subproblem for each node $n\in \nodes$, we can have one subproblem defined for each stage $t\in \horizon$ and MC state $m_t\in \markov_t$, that is, all nodes $n \in \staget$ with $m_t(n)= m_t$ will share the same cost-to-go function [Philpott  De Matos, 2012]. However, this might not be the case when we consider the aggregated variables $\za$ because
cost-to-go functions with the same MC state might be associated with different aggregated variables depending on the chosen transformation. In such a case, the cost-to-go functions are defined over two different sets of variables and, thus, the set of variables and constraints corresponding to the optimization problems are different. To see this, consider two different nodes $n,n'\in \staget$ with the same MC state (i.e., $m_t(n) = m_t(n') = m_t$) but with different sets of (aggregated) integer variables $\za$ (i.e., $\phi_t(n)\neq \phi_t(n')$), then we need to define two separate SDDP subproblems for $n$ and $n'$.
In other words, we can only define the same SDDP subproblem for two different nodes $n,n'\in \staget$ if they correspond to the same MC state (i.e., $m_t(n)=m_t(n')$) and the same transformed variables (i.e., $\phi_t(n) = \phi_t(n')$). In the following example we use a policy graph—a concept introduced by [Dowson, 2020]—that represents the relationship between SDDP subproblems. Note that the policy graphs for (<ref>) are constructed based on the MC transition graph and the chosen transformation.

Figure <ref> illustrates different policy graphs for the  , , and  transformations based on the illustrative example shown in Figure <ref>.
Recall that the underlying stochastic process is defined by an MC with two states (i.e., light 1 and dark 1). Since at a minimum we should have one SDDP subproblem defined for each MC state, we have at least two SDDP subproblems defined in each stage $t\geq 2$. In the case of , since nodes with the same MC state share the same set of aggregated integer state variables, we have exactly two SDDP subproblems per stage (i.e., one per MC state) and the policy graph only depends on the MC transitions (i.e., a Markovian policy graph). The same applies for , although we only have one set of aggregated integer state variables per stage. Since their corresponding policy graphs are identical, we use / to denote the first policy graph in Figure <ref>. On the other hand, the SDDP subproblem structure for  is quite different: nodes with the same MC state can be associated with different aggregated integer state variables, because the  transformation depends not only on the MC state in the current node but also its parent node. For example, nodes (3,111) and (3,111) are associated with the same stage and MC state ($m_{3}(\mcLight{1}\mcDark{1}\mcDark{1}) = m_3(\mcLight{1}\mcLight{1}\mcDark{1}) = \mcDark{1}$, see Figure <ref>), but since the MC states of their parent nodes differ, they have different sets of aggregated integer state variables (i.e., $\za_{3,\pentaDarkDark}$ and $\za_{3,\pentaLightDark}$, respectively). Therefore,  has one SDDP subproblem for each pair of MC states, corresponding to both the current node and the parent node for stages $t\geq 3$. The policy graph for  can be seen as an aggregated version of its scenario tree (see Figure <ref>) with one copy for each type of nodes. $\square$

SDDP policy graphs and subproblems for various types of transformations, , , and .

To formally define the SDDP subproblem for (<ref>), we represent the SDDP subproblem corresponding to node $n\in \staget$ in the scenario tree by parameterizing it with mapping $\sub(n) := (t, m_{t}(n), \phi_t(n))$. Consider $\subSett$ to be the set of all SDDP subproblems defined for stage $t\geq 2$. The number of subproblems in stage $t=2$ is always the same (i.e., $|\subSet_2| = \stage_2$)  since each node in that stage is associated to a different MC state, whereas the number of subproblems in stages $t > 2$ depends on the transformation employed. For a given SDDP subproblem $\sub\in \subSett$, we define $\subnodest(\sub)$ as the set of nodes associated with $\sub$, i.e., all nodes $n \in \staget$ such that $\sub(n) = \sub$.

To represent the relationship between SDDP subproblems illustrated by a policy graph, let $\childsub(\sub)$ be the set of children subproblems for each subproblem $\sub$ in the policy graph. Specifically, $\sub'$ is a child of $\sub$  if $\sub'\in \subSet_{t+1}$ and there exists a node $n' \in \subnodes_{t+1}(\sub')$ such that its parent is associated with $\sub$ (e.g., the children set of 1 in the  policy graph is $\childsub(\squareDark{1})=\{\pentaDarkDark, \pentaDark{1}\}$, as shown Figure <ref>). The transition probability from $\sub\in \subSett$ to one of its children $\sub' \in \childsub(\sub)$ is given by the MC transition probability, that is, $\pprob_{\sub\sub'}= \prob_{a(n')n'}$ for any node $n' \in \subnodes_{t+1}(\sub')$ with parent $a(n') \in \subnodest(\sub)$. Lastly, since there is at least one SDDP subproblem defined for every MC state, the random coefficients involved in the SDDP subproblems can be parameterized by the corresponding MC states, for example, random cost vector in $\sub\in \subSett$ is defined as $d_\sub = d_n$ for any $n \in \subnodest(\sub)$.

To preserve the stage-wise dependency structure between SDDP subproblems in two consecutive stages, we introduce local copies of the first-stage integer state variables $\za$ in the subproblems. In other words, we copy the first-stage integer solution from an SDDP subproblem to its children subproblems to ensure that each subproblem has the needed information given by $\za$. We note that the use of local copies is common in other SDDP variants to maintain the stage-wise dependence structure of the problem (see, e.g., [Zou et al., 2019]). Then, the SDDP subproblem $\subn:=\sigma(n)$ for any node $n\in \staget$ with $t>1$ is given by:

\begin{align}
\Qsub_{\subn}(\xnparentsub, \zcopyparent) = \min\; & d_{\subn}^\top \xnsub + h_{\subn}^\top \ynsub + \sum_{\subn' \in \childsub(\subn)}\pprob_{\subn \subn'} \theta_{\subn'}  \nonumber \\
\mbox{s.t.} \; & J_\sub \xnsub \geq  F_\sub \xnparentsub + f_\sub,  \\
&\zcopyn = \zcopyparent,  \label{eq:sub_zcopy} \\
& C_\sub \xnsub + D_\sub \zcopyphin + E_\sub \ynsub  \geq A_\sub \xnparentsub + B_\sub  \zcopyphiparent + b_\sub, \label{eq:sub_allvars} \\
& \theta_{\subn'} \geq \alpha_{\subn'}^\top \xnsub + \beta_{\subn'}^\top \zcopyn + \gamma_{\subn'}, \qquad \forall  (\alpha_{\subn'},\beta_{\subn'},\gamma_{\subn'})\in \coeffs_{\subn'}, \ \subn' \in \childsub(\subn), \label{eq:sddp_cut}
\end{align}

where $\zcopyn$ denotes the continuous local copy of the first-stage integer state variables $\za$, which comes from its parent subproblem $\zcopyparent$ via inequality (<ref>).  We note that since $\Qsub_{\rootnode}$ contains the original first-stage integer variables $\za$, the SDDP subproblems for stage $t = 2$ directly receive these variables instead of making local copies from the root (i.e., $\Qsub_{\sub}(\xnparentsub, \za)$ for each $\sub \in \childsub(\rootnode)$, see Figure <ref>). A variable $\theta_{\sub}$ is introduced for each SDDP subproblem $\sub$ to represent the outer approximation of the cost-to-go function, and a set $\coeffs_{\subn}$ is introduced to store the coefficients associated with all the cutting planes constructed in this approximation during the solution procedure.

Finally,  $\Qsub_{\subn}$ requires information from the parent of $\subn$ (i.e., $\parent(\sub)$), however, an SDDP subproblem can have multiple parents in the policy graph, (e.g., both  and  are parents of , see Figure <ref>). This ambiguity is not an issue for our SDDP implementation as we only consider one scenario at a time in the SDDP forward pass. So each time $\Qsub_{\subn}$ is solved, a specific node $n \in \subnodest(\sub)$ is considered as its parent node, and, thus, the parent of the SDDP subproblem would be $\parent(\subn) = \sub(\parent(n))$.

§.§.§ SDDP Algorithm.

Algorithm <ref> describes the SDDP sub-routine implemented in our B&C algorithm. The procedure takes a candidate solution $(\xrootsol,\zasol,\thetasol)$ from the master problem (<ref>) and one child node of the root $n'\in \childs(\rootnode)$, and runs the SDDP algorithm until we find a valid inequality that cuts off $(\xrootsol,\zasol,\thetasol)$ or show that no violated cuts can be found.  We use the multi-cut version of SDDP (i.e., one variable to approximate the objective value of each SDDP subproblem [Philpott  De Matos, 2012]) and  the single-scenario iteration scheme (i.e., sample one scenario at a time from the scenario tree during the SDDP forward pass [Philpott  Guan, 2008]). Lastly, our algorithm considers the quick pass version of SDDP  (i.e., we evaluate all the subproblems in a sample path before generating cuts in the backward pass [Morton, 1996]). We refer the reader to [Füllner  Rebennack, 2021] for further details on this and other SDDP variants.

1.2SDDP procedure
[1]
$(\xrootsol,\zasol,\thetasol), n', \epsilon, K,  \algExact$
Sample $K$ paths from the sub-tree rooted at $n'$ (i.e., $\tree(n')$) and store indices in $\K$.
:= .
$k \in \K$
$t \in \{2,...,\periods\}$   %% Forward pass sub-routines
Select node of scenario $k$ and stage $t$, $n:=n_{t,k}$, and its subproblem $\sub: =\sub(n)$.
Solve subproblem $\sub$, let $\Qsubsol_{\sub} :=\Qsub_{\sub}(\hat{x}_{\parent(\sub)},\zcopyparentsol)$. Save primal and dual solutions.
$t \in \{\periods,...,2\}$  %% Backward pass sub-routines
Select node of scenario $k$ and stage $t$, $n:=n_{t,k}$, and its subproblem  $\sub: =\sub(n)$.
$|\thetasol_{\sub} - \Qsubsol_{\sub}| \geq \epsilon|\Qsubsol_{\sub}|$
Add inequality (<ref>) to all subproblems with variables $\theta_{\sub}$, including $\Qsub_{\parent(\sub)}$.
:= .
$t = 2$  $(\alpha_{\sub(n')}, \beta_{\sub(n')}, \gamma_{\sub(n')})$.  %% Add cut to master problem
($K$, , ).  %% Termination criteria

We now provide a detailed explanation of the main components of the algorithm . We first sample (without replacement) $K$ scenarios (sample paths) from the sub-tree rooted at node $n'\in \stage_2$ to generate cuts for its associated variable $\theta_{n'}$ in the master problem. For each scenario, we then perform the forward and backward pass of the SDDP algorithm (lines 5-14). The forward pass iterates over all stages from $t=2$ to $t = \periods$ on a randomly sampled scenario (sample path), selects the associated node at each stage, and solves the corresponding SDDP subproblem (lines 6-8).  The backwards pass iterates over all stages from $t = \periods$ to $t = 2$ (lines 9-14) and checks the relative difference between the approximated value $\thetasol_{\sub}$ and $\Qsubsol_{\sub}$ (line 11). If this difference is larger than a given tolerance $\epsilon$, we generate a Benders cut as represented in (<ref>). In particular, the $(\alpha_{\sub'}, \beta_{\sub'}, \gamma_{\sub'})$ coefficients for any $\subn' \in \childsub(\subn)$ are given by:

\begin{align*}
\alpha_{\subn'} & = {\pi^\eqref{eq:st_onlyx}}^\top F_{\subn'} + {\pi^\eqref{eq:sub_allvars}}^\top  A_{\subn'}, \qquad \beta_{\subn'} =  {\pi^\eqref{eq:sub_allvars}}^\top B_{\subn'}, \qquad
\gamma_{\subn'}  = \Qsubsol_{\subn'} - \left( \alpha_{\subn'}^\top \hat{x}_{\subn} + \beta_{\subn'}^\top \zcopyphisol \right).
\end{align*}

where ${\pi^\eqref{eq:st_onlyx}}$ and $\pi^\eqref{eq:sub_allvars}$ are the dual solutions associated with inequalities (<ref>) and (<ref>) of subproblem $\sub'$, respectively. This cut is added to the parent subproblem (i.e., $\Qsub_{\sub(\parent(n))}$) but it is also valid for any subproblem in stage $t-1$ involving variable $\theta_{\sub(\parent(n))}$  [Philpott  De Matos, 2012]. We note that while the forward pass randomly samples scenarios rooted at one particular child node of the root node, the backward pass can add cuts to subproblems in subtrees rooted at other children nodes of the root node. Also, we consider feasibility cuts during the forward pass if one of the subproblems is infeasible.
Lastly, we terminate the current SDDP iterations for scenarios rooted at $n'$ if a cut is added to the first-stage (master) problem (i.e., $t=2$), otherwise, we continue iterating over the remaining scenarios until the termination criterion is met.

Procedure  ensures that we terminate the SDDP algorithms when we can guarantee that no cut can be added to the master problem. To do so, we increase the number of sampled scenarios $K$ to be equal to the total number of scenarios (i.e., $K=|\stage_\periods|$)) when no cut was added using the original sample size. Once $K=|\stage_\periods|$, we continue iterating until no violated cut can be added to any subproblem (i.e., =). We note that this termination criterion can lead to several iterations over the full set of scenarios of the scenario tree, which is computationally expensive. However, this is necessary to ensure an optimal solution to (<ref>).

An illustrative example of the SDDP integrated B&C algorithm. It depicts decomposition structure and the two main SDDP subroutines (i.e., forward and backward pass).

Figure <ref> illustrates the main components of our SDDP integrated B&C algorithm for the  transformation. Note that in this case an SDDP subproblem is defined for each MC state in each stage $t > 1$. The two illustrations show the decomposition structure over the policy graph: the root node 1 corresponds to the master problem while other nodes correspond to the SDDP subproblems. The drawing on the left depicts the forward pass of Algorithm <ref> for scenario $\mcLight{1} \rightarrow \mcDark{1}  \rightarrow \mcDark{1} \rightarrow \mcDark{1}$. Note that the solution of the SDDP master problem, $(\zasol, \hat{x}_{\mcLight{0.6}})$, is used as part of the input parameters to 1 and, similarly, $(\hat{\zcopy}_{\squareDark{0.6}}, \hat{x}_{\squareDark{0.6}})$ is used as part of the input to 1. Recall that the $\zcopy$ variables are just the copies of $\za$, so $\zasol=\hat{\zcopy}_{\squareDark{0.6}} = \hat{\zcopy}_{\pentaDark{0.6}}$. The drawing on the right illustrates the backward pass where the dual optimal solution of an SDDP subproblem is used to generate cuts to improve the outer-approximation of the cost-to-go function defined at its parent subproblems. For example, the optimal dual solution of 1 is used to generate cuts for subproblems 1 and 1.

§.§.§ Lower Bound Computation.

Algorithm <ref> guarantees an optimal solution to (<ref>), but the procedure can be computationally expensive because the number of SDDP subproblems depends on the underlying MC structure and the chosen transformation. In addition, we may need to call the SDDP sub-routine multiple times to correctly evaluate each candidate integer relaxation solution encountered during the B&C procedure. With this in mind, we propose an alternative to the aforementioned procedure to considerably reduce its computation time and return valuable information about the problem in the form of a lower bound. This lower bound can be used, for example, to provide an overestimate of the optimality gap associated with a known feasible solution. We utilize this lower bound to evaluate the quality of solutions obtained by our proposed 2SLDR approximations (see Section <ref> for details).

The overall idea for this lower-bounding technique is to apply a limited number of iterations of the SDDP sub-routine. First, we  increase the cut violation tolerance parameter $\epsilon$ only to add cuts that have a significant relative difference to its true value (e.g., using $\epsilon=0.1$ instead of a conservative value $\epsilon=0.0001$). Also, we limit the number of rounds in the  procedure in Algorithm <ref>. Specifically, when Boolean variable $\algExact= \algFalse$, we limit the number of rounds to only three and avoid increasing the number of sampled scenarios $K$. This is in contrast with the $\algExact= \algTrue$ version (see the previous subsection), where we iterate until no more cuts are added to any subproblem and the number of sampled scenarios increases to cover all possibilities.
These modifications to the termination criteria considerably reduce the computational time of the SDDP procedure. However, we only obtain lower-approximations to the cost-to-go functions as we might be missing cuts in the SDDP master and subproblems, resulting in a lower bound for (<ref>).

§.§ Two-stage Linear Decision Rule Approximation

In this section, we describe an approximation scheme to create feasible solutions for (<ref>) by imposing additional structures to the decision policies in the form of 2SLDR  [Bodur  Luedtke, 2022]. Imposing additional structures in such a way transforms the multi-stage stochastic program into a two-stage stochastic program, which requires significantly less computational effort than the B&C framework integrated with the SDDP algorithm that we presented in the previous section. Moreover, the 2SLDR can be customized to leverage the structure of the underlying stochastic process, and we provide novel schemes for the MC case.

Linear decision rules (LDRs) have been widely used in optimization under uncertainty [Garstka  Wets, 1974], most notably in adjustable robust optimization [Ben-Tal et al., 2004]. LDRs have also been employed in multi-stage stochastic linear programs (MSLPs) to compute primal and dual bounds [Shapiro  Nemirovski, 2005, Kuhn et al., 2011] by transforming the multi-stage problem into a single-stage problem. Recently, [Bodur  Luedtke, 2022] propose an LDR framework to transform an MSLP into a two-stage stochastic linear program by applying LDRs to only the state variables—a procedure we refer to as 2SLDR—and is guaranteed to produce superior decision policies than standard LDRs.

We now provide details on how to apply 2SLDRs to (<ref>). Let $\xi_t$ be the vector of random variables at stage $t\in \horizon$, and $\xi_{t,n}$ be the realization of these random variables at node $n\in \staget$,  i.e.,
$ \xi_{t,n} = \{A_n,B_n,C_n,D_n,E_n,F_n,J_n,b_n,c_n,d_n,f_n,h_n\}. $ Note that $\xi_\rootnode$, the data associated with the root node of the scenario tree, is assumed to be a deterministic vector. In addition,  $\xi^t=\{\xi_1,...,\xi_t\}$ and $\xi^t_n=  \{\xi_{1,\rootnode},...,\xi_{t-1, \parent(n)}, \xi_{t,n}\}$ represent the trajectory of the random variables up to stage $t$ and their realizations for a particular node $n$ in the scenario tree, respectively. We consider $\xi_t$ and $\xi^t$ as vectors of random variables with dimensions $\dimxit$ and $ \dimxitt =\sum_{i = 1}^t \dimxit$ for all $t\in \horizon$, respectively, and the same applies to the realizations of these random variables.

The general idea behind 2SLDRs is to replace the original state variables $\xn$ with a new set of variables for each node $n \in \nodes$ by means of a linear relation of node-specific data. By doing so, we impose additional (linear relation) structure over the state variables, and in return we transform the multi-stage stochastic program into a two-stage stochastic program [Bodur  Luedtke, 2022]. Specifically, 2SLDRs employ basis functions $\mapldr(\cdot)$ to map $\xi^t_n$ into values that are then used to create a linear representation for the state variables $\xn$, i.e., $\xn = \varldr_t^\top \mapldr(\xi^t_n) \in \R^k$, where $\varldr_t$ is a vector of the so-called LDR variables at stage $t$. Although the basis functions $\mapldr(\cdot)$ can be any form of functions that map the node-specific data $\xi^t_n$ to certain values, for ease of exposition, we consider simple basis functions that return the actual random variable realizations, that is, $\mapldr(\xi^t_n) = \xi^t_n$. In this case, the LDR is given by:

\begin{equation}\label{eq:ldr}
\xn = \varldr_t^\top \xi^t_n, \ \forall n \in \staget, t \in \horizon.
\end{equation}

We observe that these simple basis functions return high-quality decision policies for our case study (see Section <ref>), although more complex basis functions can potentially lead to better decision policies [Chen et al., 2008, Bampou  Kuhn, 2011].

§.§.§ Aggregation Framework for 2SLDR.

Similar to the aggregation framework introduced in Section <ref>, we create different aggregation  schemes for the LDR variables $\varldr$, and as a result, obtain approximate 2SLDR-based decision policies. Table <ref> summaries a few LDR variable structures that we consider. The first two LDR variants ( and ) are common in the literature [Shapiro  Nemirovski, 2005, Kuhn et al., 2011], while  is a novel LDR variant based on the MC structure of the stochastic process.

Name     LDR variables for $t\in \horizon$     Standard basis function     Resulting LDR

$\varldr_t \in \R^{k\times \dimxitt}$     $\mapldr(\xi^t_n) = \xi^t_n$     $\xn =  \varldr_t^\top \xi^t_n$

$\varldr_{t} \in \R^{k\times \dimxit}$     $\mapldr(\xi^t_n) = \xi_{t,n}$     $\xn = \varldr_t^\top \xi_{t,n}$

$\varldr_{t,m}\in \R^{k\times\dimxit},  \quad \forall\; m \in \markov_t$     $\mapldr(\xi^t_n) = \xi_{t,n}$      $\xn =  \varldr_{t, m(n)}^\top \xi_{t,n}$

Examples of LDR variable structures and basis functions for each $n\in \staget$ and $t\in \horizon$.

The first LDR variant, , which we refer to as a stage-history LDR, constructs the linear relation based on the trajectory of realizations of the random variables up to the current stage. The  variant can lead to high-quality decision policies (thanks to the amount of information incorporated) but may result in a large number of LDR variables, especially for problems with a large number of stages. An alternative is  which constructs the linear relation only based on the realization of the random vector at the current stage, which leads to a smaller number of LDR variables at the expense of lower solution quality [Bodur  Luedtke, 2022].

In addition to these LDR variants, we  propose  which takes advantage of the MC structure of the underlying stochastic process. The idea is to have one copy of LDR variables $\varldr_{t,m}$ for each MC state $t$ and stage $m$. It is clear that  will lead to better policies compared to  because both LDR variants use the same information (realizations of random vectors), but the latter has a single copy of the LDR variables for all nodes in the same stage. We note that, however,  and  are incomparable in general, as shown in our empirical results.

§.§.§ 2SLDR Models.

We now present the resulting model after applying the LDR transformation (<ref>) to (<ref>). For ease of exposition, we only present the 2SLDR variant  and the other two alternatives are similar. The overall idea is to replace the state variables $\xn$ with the LDR variables using linear function (<ref>). The resulting model for the root node is as follows:

\begin{align}
\Qldr_{\rootnode} = \min\; & %\left\{(\yroot, \varldr, \za) \in \constldr_{\rootnode}: \;
\sum_{n\in \nodes}\prob_n \left( c_n^\top \zphi + d_n^\top(\varldr_{t(n)}^\top\xi^{t(n)}_n) \right) + h_\rootnode^\top \yroot + \sum_{n \in \nodesnoroot }\prob_{n} \Qldr_n(\varldr,\za) \tag{$P^L$}  \label{mod:ldr} \\
\text{s.t.}\; &  H_{\rootnode} \zroot \geq g_{\rootnode}, \;
J_{\rootnode} \varldr_{1}^\top\xi^{1}_\rootnode \geq f_{\rootnode}, \;
C_{\rootnode} \varldr_{1}^\top\xi^{1}_\rootnode + D_{\rootnode} \zroot + E_{\rootnode} \yroot  \geq  b_{\rootnode}, \label{eq:ldr-constr0} \\
& H_n \zphi \geq  G_n \zphiparent + g_n, & \forall n \in \nodesnoroot,  \label{eq:ldr-constr1} \\
& J_n \varldr_{t(n)}^\top\xi^{t(n)}_n \geq F_n  \varldr_{t(a(n)) }^\top\xi^{t(a(n))}_{\parent(n)}, &  \forall n \in \nodesnoroot, \label{eq:ldr-constr2} \\
& (\yroot, \varldr, \za)  \in  \R^r \times \R^{k\cdot \sum_{t \in \horizon}\dimxitt} \times \Z^{\ell \cdot \sum_{t \in \horizon} q_t}. \nonumber
\end{align}

Model (<ref>) contains all the root node variables and constraints of (<ref>) and includes all the LDR variables and its associated constraints (<ref>) from all nodes in the scenario tree. We note that there is no need to define LDR variables for the root node since the $\xroot$ variables are already first-stage variables, but we include them in the model for ease of exposition. Lastly, the cost-to-go function for node $n \in \nodesnoroot$ only depends on the first-stage variables $\varldr$ and $\za$ and is given by:

\begin{equation*}
\Qldr_n(\varldr, \za) = \min_{\yn \in \R^r} \left\{ h_n^\top \yn \mid  C_n \varldr_{t(n)}^\top\xi^{t(n)}_n + D_n \zphi + E_n \yn  \geq A_n \varldr_{t(n) -1 }^\top\xi^{t(n)-1}_{\parent(n) } + B_n  \zphiparent + b_n\right\}.
\end{equation*}

Note that model (<ref>) is indeed a two-stage stochastic programming model. By applying (<ref>) to the state variables, we remove the dependency between parent and child nodes and, thus, all cost-to-go functions are only functions of the first-stage decision variables. This two-stage stochastic program is easier to solve than its multi-stage stochastic programming counterpart since it eliminates the nested dependency between the cost-to-go functions in different stages.

Formulation (<ref>) considers one cost-to-go function for each node in the scenario tree. However, the number of cost-to-go functions can be considerably smaller depending on the LDR variant and transformation $\Phi_t$. For example,  requires one cost-to-go function for each node in the scenario tree because it considers the entire trajectory of realizations of random variables. Therefore, the cost-to-go functions for  are unique for each node even if they share the same stage and MC state. On the contrary,  and  only employ local information, so the cost-to-go functions for two nodes $n, n'\in \staget$ are the same if $\sub(n) = \sub(n')$, that is, if they share the same stage $t$, MC state $m(n)=m(n')$, and aggregated variables $\phi_t(n)= \phi_t(n')$. Therefore, we have the same number of distinct cost-to-go functions as SDDP subproblems, which can be significantly smaller than the size of the scenario tree. In conclusion,  and  have a significant computational advantage over  due to the fewer LDRs variables and cost-to-go functions employed (along with all the variables and constraints necessary to define them).

We solve the resulting two-stage stochastic programming model using Benders decomposition, a standard approach in the stochastic programming literature (see, e.g., [Zverovich et al., 2012]). Our specific implementation details can be found in Appendix <ref>.

§ CASE STUDY: A HURRICANE DISASTER RELIEF PLANNING PROBLEM

In this section, we present a case study to exemplify the aggregated framework and evaluate the proposed solution methodologies. We consider a novel class of multi-period hurricane disaster relief logistics planning (HDR) problem, where a number of contingency modality options can be activated during the planning horizon to increase the capacities at the distribution centers (DCs).

The HDR problems are well-studied in the operations research literature because of their practical impacts on reducing economic loss and human suffering [Graumann et al., 2006, Seraphin, 2019]. The operational challenges in HDR problems arise from the hurricanes' stochastic nature [Sabbaghtorkan et al., 2020]. To address these challenges, most of the existing literature on HDR problems focuses on two-stage stochastic programming models, where the first-stage model corresponds to strategic facility location and resource pre-positioning decisions in terms of shelter preparation and resource allocation, while the second-stage model deals with disaster relief logistics decisions once the damage and demand information is revealed after the hurricane's landfall (see, e.g., [Duran et al., 2011, Lodree Jr et al., 2012, Davis et al., 2013, Alem et al., 2016], and a survey paper by [Sabbaghtorkan et al., 2020]).

Recent works have also considered the multi-stage HDR variants where sequential logistics decisions are made in a rolling horizon (RH) fashion [Pacheco  Batta, 2016, Siddig  Song, 2022, Yang et al., 2022]. [Pacheco  Batta, 2016] consider a relief supply preposition and reposition problem starting from the time when the hurricane is first detected, with updated forecast information on the hurricane's attributes every six hours. The authors proposed a forecast-driven dynamic model for the problem, and a solution approach combining scenario analysis and RH. Similarly, [Siddig  Song, 2022] consider a multi-period preposition problem where the demand is realized at landfall (i.e., last period). The authors proposed an MSLP formulation and showed the value of multi-stage stochastic programming solutions compared to solutions given by an RH approach and a two-stage approximation. Lastly, [Yang et al., 2022] focus on the diesel fuel supply chain problem for hurricane relief logistics, applying a two-stage stochastic programming model with updated forecast information at each roll of an RH procedure.

Similar to [Pacheco  Batta, 2016] and [Siddig  Song, 2022], the class of HDR problems that we study is concerned with distributing relief commodities to the shelters over multiple stages. The main difference in our HDR problem is that we consider a set of contingency modality options that can be activated during the planning horizon to increase the inventory capacities of the DCs. In addition, we consider that the demand realization starts at the beginning of the planning horizon, which is assumed to be the time when the hurricane starts to form at sea, while previous works consider demand realization only at the hurricane's landfall.

Despite their wide range of problem settings, the HDR problems have several assumptions in common. First, the hurricane movement and intensity are commonly modeled using an MC  [Taskin  Lodree Jr, 2010, Pacheco  Batta, 2016] and the estimated demand at each shelter can be modeled using a deterministic mapping from the current hurricane state (i.e., hurricane location and intensity)[Siddig  Song, 2022]. Once the hurricane is first observed, it is common to consider that each decision epoch corresponds to a six-hour time window, which is consistent with the frequency of the hurricane forecast update by the National Hurricane Center [Pacheco  Batta, 2016, Regnier  MacKenzie, 2019]. The HDR problem we study follows these assumptions as well, and we next introduce more details about the problem.

§.§ Problem Description and Model Formulation

Our novel multi-stage HDR problem is concerned with the logistics decisions of producing and distributing relief commodities to shelters prior to the landfall of a hurricane event. These relief commodities are used to fulfill the demand for civilians that are evacuated to shelters. As in related HDR problems, we consider that the demand can be estimated according to the location and intensity of the hurricane, which are updated every six hours. We assume that shelters, DCs, and possible positions of the hurricane are located in a two-dimensional grid with a discrete set of $x$- and $y$-coordinates. The main goal is to produce and transport relief commodities from the DCs to the shelters in order to minimize the total cost, which consists of the penalty cost for unsatisfied demand and operational costs incurred in all time periods.

In addition, we consider a set of contingency modalities that can be activated to increase the DCs' capacities. These modalities make the disaster relief logistics system more resilient to extreme hurricane events that may lead to unexpectedly high demand at the price of higher logistics costs. We make a practically relevant assumption that these modalities can be activated at most once during the planning horizon, and once activated, the modality will stay the same until the end of the planning horizon. Therefore, an active modality represents an incremental capacity increase at every stage on the planning horizon. These considerations are consistent with the practical constraints associated with these large-scale critical logistics operations (usually military aids are involved), which have this “all or nothing” feature and have to be implemented in small increments. Lastly, we consider a one-leg delay for these modalities (i.e., the capacity starts to increase from the stage immediately after the stage when a modality activation decision is made), which partially captures the logistical challenges of implementing changes during the disaster relief operations.

We assume that the hurricane's evolution can be characterized by an MC $\markov$. Each MC state corresponds to two attributes of the hurricane: location and intensity, that is, $m=(m^x, m^y, m^i)\in \markov$ where $m^x$ and $m^y$ are the $x$- and $y$-coordinates of the hurricane's location, respectively, and $m^i$ represents the hurricane's intensity level. We consider independent probability transition matrices for each hurricane attribute. We employ the transition probability matrix described in [Pacheco  Batta, 2016] for the hurricane intensity. For the hurricane movement, we assume that the hurricane originates from the bottom row of the grid (i.e., $m^y = 0$ for the initial MC state) and that $m^y$ increases by one in each period (i.e., the hurricane advances upwards by one step in each period) and will reach land in exactly $T$ periods. See Appendix <ref> for further details.

Given the stochastic behavior of the hurricane, we consider an MSILP model with a stage defined for each time period when new information is realized (i.e., the hurricane moves to a new MC state). We construct a scenario tree $\tree$ by considering  an initial MC state and all possible MC transitions for up to $T$ stages (i.e., until  landfall). We consider a set of contingency modalities $\modality$, a set of shelters $\shelters$, and a set of DCs $\centers$, where each DC $j\in \centers$ has an initial capacity $C_j$ and inventory $I_j$. At each node $n\in \nodes$ of the scenario tree $\tree$, one needs to decide the amount of relief commodities $v_{j}\geq 0$ to produce for each DC $j \in \centers$, the amount of relief commodities $y_{ij}\geq 0$ to transport from each DC $j\in \centers$ to each shelter $i \in \shelters$, and the amount of unsatisfied demand $w_i \geq 0$ at each shelter $i\in \shelters$. We introduce two continuous state variables $x_n = (\xinv_{n}, \xcap_{n})$ to represent the inventory and capacity for each DC $j\in \centers$ at node $n\in \nodes$, respectively. Lastly, we introduce a binary state variable $z_{n\ell} \in \{0,1\}$ to represent whether or not a contingency modality $\ell \in \modality$ is active at node $n \in \nodes$.

Each node $n\in \nodes$ is associated with realizations of random variables according to the corresponding MC state of the hurricane. Specifically, for a given node $n \in \staget$, $t \in \horizon$, the demand of a shelter $i\in \shelters$, $d_{ni}$, the production cost at DC $j\in \centers$, $q_{nj}$, and the transportation cost from DC $j$ to shelter $i$, $f_{nij}$, depend solely on the MC state $m_t(n)$ associated to node $n$. Thus, two nodes $n,n'\in \staget$ with the same MC state at stage $t$, $m_t(n)=m_t(n')$, share the same realizations of random variables. The remaining model parameters are assumed to be deterministic, which include the unit inventory cost $g_j$ at $j\in \centers$, the contingency modality cost $c_\ell$ for $\ell\in \modality$, the penalty cost for each unit of unsatisfied demand $b_i$ for $i\in \shelters$, and the capacity increase $K_{j\ell}$ at DC $j$ in each stage when modality $\ell$ is active. Further details on the problem instance  generation can be found in Appendix <ref>. The MSILP model for the HDR problem is given by:

\begin{align}
\min \quad &  \sum_{j \in \centers}\left( g_{j}\xinv_{\rootnode j} + q_{\rootnode j} v_{j} +\sum_{i \in \shelters}  f_{\rootnode ij} y_{ij} \right) + \sum_{i \in \shelters} b_{i}w_{i} + \sum_{\ell \in \modality} z_{\rootnode \ell}c_{\ell} + \sum_{n \in \childs(\rootnode)} \lefteqn{\pprob_{\rootnode n} Q_n(\xroot, \zroot) } \tag{$\textit{HDR}$} \label{dr_root} \\
\text{s.t.} \quad & \sum_{j \in \centers} y_{ij} + w_{i} \geq d_{\rootnode i}, & \forall i \in \shelters,  \label{dr_root:demand} \\
& \xinv_{\rootnode j} = I_{j} - \sum_{i \in \shelters} y_{ij} + v_{j}, & \forall j \in \centers, \label{dr_root:inventory} \\
& v_{j} \leq \xcap_{\rootnode j},  & \forall j \in \centers, \label{dr_root:production}\\
& \xcap_{\rootnode j} = C_j,  & \forall j \in \centers, \label{dr_root:capacity} \\
& \sum_{\ell \in \modality} z_{\rootnode \ell} \leq 1, \label{dr_root:modality} \\
& \xinv_{\rootnode j}, \xcap_{\rootnode j},   v_{j}, w_{i}, y_{ij} \geq 0, \; \zroot{}_\ell \in \{0,1\}, & \forall j \in \centers, \; i \in \shelters, \; \ell \in \modality.  \nonumber
\end{align}

Constraints (<ref>) and (<ref>) represent the demand and inventory constraints, respectively. Constraints (<ref>) enforce the production capacity and constraints (<ref>) represent the initial capacity of each DC. Constraint (<ref>) enforces that at most one modality can be activated at the root node of the scenario tree. Lastly, the cost-to-go function for a node $n\in \nodesnoroot$ is given by:

\begin{align}
Q_n(\xnparent, \znparent) &  =  \min \sum_{j \in \centers}\left( g_{j}\xinv_{n j} + q_{n j} v_{j} + \sum_{i \in \shelters}f_{nij} y_{ij} \right) +  \sum_{i \in \shelters} b_{i}w_{i} + \lefteqn{\sum_{\ell \in \modality} z_{n \ell}c_{\ell}  %\nonumber \\
+\!\! \sum_{n' \in \childs(n)}\!\!\pprob_{n n'} Q_{n'}(\xn, \zn) } \quad
\nonumber \\
\text{s.t.} \quad  & \sum_{j \in \centers} y_{ij} + w_{i} \geq d_{ni}, & \forall i \in \shelters, \label{dr_node:demand} \\
& \xinv_{nj} = \xinv_{a(n)j} - \sum_{i \in \shelters} y_{ij} + v_{j}, & \forall j \in \centers, \label{dr_node:inventory} \\
& v_{j} \leq \xcap_{nj},  & \forall j \in \centers,\label{dr_node:production} \\
& \xcap_{nj} = \xcap_{a(n)j} + \sum_{\ell \in \modality} K_{j \ell}z_{a(n) \ell},  & \forall j \in \centers, \label{dr_node:capacity}\\
& \sum_{\ell \in \modality} z_{n\ell} \leq 1, &  \label{dr_node:modality}\\
& z_{a(n)\ell} \leq  z_{n\ell}, & \forall \ell \in \modality,  \label{dr_node:modality2}\\
& \xinv_{n j}, \xcap_{n j},   v_{j}, w_{i}, y_{ij} \geq 0, \; \zn{}_\ell \in \{0,1\}, & \forall j \in \centers, i \in \shelters, \ell \in \modality. \nonumber
\end{align}

The meaning of constraints (<ref>)-(<ref>) and (<ref>) are identical to their corresponding constraints at the root node. Constraints (<ref>) represent the capacity increase at a DC depending on whether or not a contingency modality is active in the previous stage. Constraints (<ref>) enforce that an active contingency modality in the previous stage will remain active in the current stage.

We note that model (<ref>) can be used for other contexts besides hurricane disaster relief. The model can represent any type of multi-stage logistics planning problem with any contingency modality structure, after a minor modification to the model formulation. However, in order to directly apply the proposed aggregated framework described in Section <ref>, the underlying stochastic process must be represented by an MC.

§.§ Aggregation Framework

We now present how to transform the MSILP model (<ref>) with integer state variables at each node of the scenario tree to an MSILP model with integer variables only at the root node, following the ideas presented in Section <ref> and Section <ref>. The general idea is to consider all integer state variables as first-stage variables and apply a suitable transformation $\maptmarkov$ to aggregate the number of variables to a manageable size. In what follows, we present the resulting aggregated model (<ref>) for a generic transformation and then show four alternatives given the specific structure of (<ref>).

Following the notation introduced in Section <ref>, we consider $\za$ to be the set of aggregated integer state variables for each node. There are two aspects to consider when applying this transformation. First, constraints (<ref>) and (<ref>) are now first-stage constraints. Second, we can omit state variables $\xcap_n$  because these variables are only used to model the capacity dependency from one stage to the next, which is fully dictated by the first-stage variables $\za$. The resulting aggregated model is:

\begin{align}
\min \ &  \sum_{j \in \centers}\left( g_{j}\xinv_{\rootnode j} + q_{\rootnode j} v_{j} + \sum_{i \in \shelters}f_{\rootnode i j} y_{ij} \right) + \sum_{i \in \shelters} b_{i}w_{i} + \sum_{n \in \nodes}\sum_{\ell \in \modality} \lefteqn{
\za_{\phi_t(n) \ell} c_{\ell} + \sum_{n \in \childs(\rootnode)}\pprob_{\rootnode n} \Qaggre_n(\xroot^I, \za) } \tag{$\textit{HDR}^A$} \label{dr_aggre} \\
\text{s.t.} \ & \eqref{dr_root:demand}-\eqref{dr_root:inventory}  \nonumber \\
& v_{j} \leq C_j,  & \forall j \in \centers, \label{dr_aggre:capacity}\\
& \sum_{\ell \in \modality} \za_{\phi_t(n) \ell} \leq 1, & \forall n\in \nodes, \label{dr_aggre:modality1} \\
& \za_{\phi_{t-1}(a(n)) \ell} \leq  \za_{\phi_t(n) \ell}, & \forall \ell \in \modality, n \in \nodesnoroot, \label{dr_aggre:modality2} \\
& \xinv_{\rootnode j}, v_{j}, w_{i}, y_{ij} \geq 0, \; \za_{\phi_t(n) \ell} \in \{0,1\}, & \forall j \in \centers, \; i \in \shelters, \ell \in \modality, n \in \nodesnoroot, \nonumber
\end{align}

where the cost-to-go function for a node $n\in \nodesnoroot$ is given by:

\begin{align}
\Qaggre_n(\xnparent^I, \za)  =  \min \ & \sum_{j \in \centers}\left( g_{j}\xinv_{n j} + q_{n j} v_{j} + \sum_{i \in \shelters} f_{nij}y_{ij} \right) +  \sum_{i \in \shelters} b_{i}w_{i} + \lefteqn{\sum_{n' \in \childs(n)}\pprob_{n n'} \Qaggre_{n'}(\xn^I, \za)} \nonumber \\
\text{s.t.}\; & \eqref{dr_node:demand}-\eqref{dr_node:inventory}  \nonumber \\
& v_{j} \leq C_j + \sum_{n' \in \patht(n)}\sum_{\ell \in \modality} K_{j \ell}\za_{\phi_t(n') \ell},  & \forall j \in \centers, \label{dr_aggre:capacity2}\\
& \xinv_{n j}, v_{j}, w_{i}, y_{ij} \geq 0, & \forall j \in \centers, \; i \in \shelters. \nonumber
\end{align}

Note that the right-hand side of constraint (<ref>) represents the current capacity of a DC by considering the modality activation and capacity increases at previous stages. Also, although constraints (<ref>) and (<ref>) are imposed for each node, the fact is that some of them might be identical because they share the same problem parameters and aggregated variables. For example, transformation  only needs one set of constraints (<ref>) per stage because of the stage-based variable aggregation.

We consider the four MC-based transformations introduced in Table <ref>, that is, , , , and , for the aggregated model (<ref>). As described in Section <ref>, transformation  considers the same set of contingency modality decisions for all nodes in a given stage; transformation  aggregates variables based on the MC state information in a given stage, which corresponds to the location and intensity of the hurricane; and transformation  incorporates the MC states at the current node and its parent node. The proposed transformation  can be customized to different settings depending on how one defines the partial MC state information. In our experiments, we consider the full MC state information of the hurricane (i.e., location and intensity) at the current node and only the hurricane intensity information at its parent node. In other words, the main difference between  and  for the HDR problem is that the former incorporates the hurricane location information at the parent node into the aggregation, while the latter does not.

§.§ Two-stage LDR approximation

As explained in Section <ref>, we can apply a further restriction to the aggregated MSILP model (<ref>) to obtain a two-stage model by applying LDR to the continuous state variables (i.e., inventory variables $\xinv_n$). We present three LDRs variants similar to the ones introduced in Table <ref>.

First, we introduce notation $\xi_{t,n}=(d_n,q_n,f_n)$ to represent the vector of realizations of random variables for each node $n \in \staget$ in the scenario tree, including the demand $d_n$, production cost $q_n$, and transportation cost $f_n$. Since the demand is one of the most important parameters for the HDR problem, in our experiment we choose to use basis functions (for constructing LDRs) that are only based on either the demand realization at the current node $\mapldr(\xi_{n}^t) = d_n$ or the entire history of demand realizations up to the current node $\mapldr(\xi_{n}^t) = (d_\rootnode,...,d_{a(n)}, d_n)$. Then, the LDR variants for the inventory variables for a given node $n\in \staget$ and DC $j \in \centers$ are defined as:

\[
\thldr: \; \xinv_{n j} = \!\!\sum_{n' \in \patht(n)}\sum_{i \in \shelters}\varldr_{t(n')tji}d_{ni} \quad\;
\tldr: \;  \xinv_{n j} = \sum_{i \in \shelters}\varldr_{t(n)ji}d_{n'i} \quad\;
\mldr: \; \xinv_{n j} = \sum_{i \in \shelters}\varldr_{t m_t(n) ji}d_{ni}
\]

where $\varldr_{t't} \in \R^{|\centers|\times |\shelters|}$ for $t' \in [t]$, $\varldr_{t} \in \R^{|\centers|\times |\shelters|}$, and $\varldr_{tm}\in \R^{|\centers|\times|\shelters|}$  for $m \in \markov_t$ are the LDR variables associated to each variant, respectively. Note that  considers the full demand history, while  and  only consider the demand realization at the current node. Also,  has one set of variables for each MC state (i.e., hurricane location and intensity), while  and   employ one set of variables per stage. We discuss the implications of these differences in Section <ref> and also in our numerical result in Section <ref>.

These LDR variants lead to three different two-stage approximations of (<ref>), which are solved via Benders decomposition as discussed in Section <ref> and Appendix <ref>. Note that while the MSILP problem has the relatively complete recourse property, these two-stage approximations do not necessarily have this property. For instance, LDR variables $\varldr$ can take negative values which may lead to an infeasible solution with negative inventory values. Thus, we incorporate Benders feasibility cuts for this and other cases in our implementation.

§ NUMERICAL EXPERIMENT RESULTS

We now present our numerical experiment results for the aggregation framework proposed in Section <ref> and solution approaches proposed in Section <ref>. We test the proposed models and approaches over a wide range of instances of the HDR problem described in Section <ref>. In what follows, we first provide details on the test instances and the experimental setup. We then present experimental results comparing the performances of different transformations for our aggregation framework in terms of the solution quality of their corresponding decision policies and their computational efforts. Finally, we evaluate the computational performances of the proposed exact method (i.e., the B&C procedure integrated with SDDP) and approximation methods (i.e., 2SLDR solutions with SDDP-based lower bounds).

§.§ Experimental Setup

We generate instances for our  MSILP HDR problem (see Section <ref>) for two different grid sizes, $4\times5$ and $5\times6$, which correspond to mid-size and large-size instances, respectively. We experiment with different levels of initial capacities for the DCs (i.e., $C_j$ for $j \in \centers$) to analyze the impact of the aggregation framework for systems with different capabilities to address demand fluctuations—DCs with higher capacities are more capable of reacting to fluctuating demand. Specifically, we consider instances where the initial capacities across all DCs are 20%, 25%, or 30% of the maximum demand across all shelters (the maximum demand is a predefined parameter for instance generation, see Appendix <ref> for additional details). We also consider two different types of contingency modality options: a conservative one (Type-1) and an aggressive one (Type-2). For instances with Type-1 modality, we consider DC capacity expansion options to be 10%, 20%, 30% or 40% per stage, and for instances with Type-2 modality, we consider DC capacity expansion options to be 15%, 30%, 45% or 60% per stage. In all, our experiments evaluate six different instance configurations for each grid size. We randomly generate a number of instances for each instance configuration, from which we randomly pick ten where the optimal decision policies corresponding to transformations  and  are different from each other. This way, we ensure that different transformations' advantages and disadvantages are reflected in their corresponding solutions in these instances.

The code for the decomposition algorithms was implemented in  using solver IBM ILOG CPLEX 20.1 with callback functions. All experiments were run using a single core/thread and a six-hour time limit. The experiments were run in the University of Toronto SciNet server Niagara, using cores with 4GM RAM each[See <https://docs.scinet.utoronto.ca/index.php/Niagara_Quickstart> for further server specifications.]. We will make our code and instances available upon publication.

§.§ Solution Quality and Computational Effort for Different Transformations

We first compare the performances of models resulting from the four transformations proposed in our aggregation framework (, , , and ) and the original, non-aggregated model (). All models are solved using the extensive form (i.e., no decomposition is used) to have a fair comparison in terms of the computational effort. These sets of experiments show results for the mid-size instances (i.e., $4\times5$ grids) because large-size instances cannot even be loaded with CPLEX due to the size of the extensive model and limited memory. We also include a sensitivity analysis discussion regarding different instance configurations (i.e., modality types and initial capacity levels).

Table <ref> summarizes the performances of the considered models in terms of the resulting objective values on instances with two types of contingency modalities and three levels of initial capacities. Columns “% Gain over ” present the
percentage of the gap between  and  that is closed by each transformation, averaged over all the instances, where the gap closed for each instance is calculated in terms of the objective value ($Obj$) as $(Obj_{\hn} - Obj_{i})/(Obj_{\hn}- Obj_{\fh}),  \forall i \in \{\ma, \pa, \mm\}$.

Policy quality comparison for all transformations in our aggregation framework.

5c|Average Objective Value     3c% Gain over

Modality     Cap.      1c     1c     1c     1c     1c|

3[1]*Type-1     20%        104,162        102,572        92,924        92,924        82,193      7.5       51.9      51.9

25%          73,879          73,442        66,650        66,238        62,972      3.2       69.4      73.1

30%          48,970          48,951        47,425        47,314        47,117      0.3       89.5      96.7

3[1]*Type-2     20%        104,135        102,625        92,612        92,612        81,182      6.7       50.9      50.9

25%          73,919          73,487        66,573        65,958        63,654      3.3       75.6      81.2

30%          48,970          48,968        47,540        47,349        47,253      0.0       80.3      97.0

First, we focus on the difference between  and  across all test instances. We see that the difference between the  and  objective values is smaller when the initial DC capacities are higher—about $4\%$ when the initial capacity is 30%, and about $28\%$ when the initial capacity is $20\%$. This observation suggests that even the most restrictive transformation  can lead to high-quality policies when the DCs are equipped with larger capacities and are more capable of addressing random demand. This is, however, not the case when the DCs have smaller capacities, where less restrictive transformations are needed for more flexibility to address random demand better. We also note that the modality type

has very limited influence on the performance of the policies.

Comparing the performances of different transformations, we see that the inclusion of information from the previous stage (i.e.,  and ) leads to a significant improvement—closing more than 50% of the gap between  and  in all cases. In particular, the objective values of these two transformations are very close to  those obtained by  when the initial capacity is 30%, but this is not the case when the initial capacity is 20%. The performances of  and  are similar, suggesting that the intensity information of the previous stage captures most of the valuable information needed. We also observe that the performances of  and  are similar (with $< 10\%$ difference) in all instances, indicating that using only the information in the current stage is insufficient to obtain high-quality policies for our test instances.

Figure <ref> illustrates the computational effort to obtain an optimal solution for each of these transformations. The plot shows the accumulated number of instances solved to optimality over the six-hour time limit (in logarithmic scale). We can see that the aggregated models require  significantly less computational effort because they solve
all 60 instances in at most an hour and a half, while   only solves 27 out of 60 instances to optimality within the time limit.

[ legend style=at=(1.06,1),anchor=north, font=, legend cell align=left,
xlabel=Time in seconds (logarithmic scale),
ylabel= # Solved,

ytick=0,  20,  40,  60,

[color=orange, mark=square, mark options=scale=0.8] table[x=HN, y=ycoord] figures/plots/time_extensive.txt;

[color=violet!70!white, mark=o, mark options=scale=0.8] table[x=MA, y=ycoord] figures/plots/time_extensive.txt;

[color=blue, mark=diamond, mark options=scale=1] table[x=MI, y=ycoord] figures/plots/time_extensive.txt;

[color=cyan, mark=x, mark options=scale=1] table[x=MM, y=ycoord] figures/plots/time_extensive.txt;

[color=green!60!black, mark=triangle, mark options=scale=1] table[x=WA, y=ycoord] figures/plots/time_extensive_WA.txt;

Number of instances solved over time. Extensive model, medium-size instances ($4\times5$ grid).

In sum, the computational results presented in this subsection show a clear trade-off between policy quality and computational effort: more restrictive transformations require less computational effort but lead to poorer quality policies. For our HDR test instances, transformation —which takes moderate information from the hurricane history—appears to have the best trade-off as it obtains high-quality policies with reasonable computational effort.

§.§ Results for Exact Methods

We now compare two solution methods for the aggregated models with four types of transformations, that is, solving the extensive form model directly with CPLEX () versus the B&C procedure integrated with the SDDP algorithm proposed in Section <ref> (). Table <ref> shows the average solution time (for instances where optimal solutions are obtained within the time limit) and the average optimality gaps (for unsolved instances) for the mid-size instances. We note that  solves to optimality all the instances for all transformations, while  fails to solve one instance for  and solves only one instance to optimality for .

From Table <ref> we see a clear advantage of  over  in both computational time and number of instances solved within the time limit across all transformations and instance configurations. These results are mostly explained by the fact that approximately 90% of the computational time is spent inside the SDDP sub-routine. One main factor is the number of SDDP subproblems for each transformation, e.g.,  and  have on average $69$ SDDP subproblems, while  has $144$ and  has $294$. This explains the relatively poor performance of .

Performance comparison between  and  for mid-size instances ($4\times5$ grid size).

8c|Average Time (sec)                  2cGap (%)

2c|     2c|     2c|     2c|     1c     1c

Modality     1l|Cap.     1c     1c|     1c     1c|     1c     1c|     1c     1c|     1c     1c

3[2]*Type-1     1l|20%              260           1,221                712         5,969             961         4,553        1,853      -        -         14.9

1l|25%              132              195                361         1,313             487         2,760           626      -        -         3.5

1l|30%                73              102                109            272             222         3,380           367      -        -         2.1

3[2]*Type-2     1l|20%              253           1,138                919         5,880          2,047       13,224        2,951       20,699      19.9      17.8

1l|25%              167              140                393         1,324             651         3,635           817      -        -         6.2

1l|30%                82              102                119            205             271         4,347           556      -        -         2.3

2c|Average              161              483                436         2,494             773         5,316        1,195       20,699         19.9        7.8

The above results may suggest that  requires less computational effort for mid-size instances than . However, as previously mentioned, the large-size instances cannot even be loaded into the solver due to the large number of variables and constraints in the formulation, so  is our only resort for these instances. Table <ref> shows the number of instances solved to optimality, the number of instances where the algorithm found a feasible solution (but could not prove optimality), and the average optimality gap over instances with a feasible solution for  over large-size instances. We see that  performs quite well for transformations  and   solving several instances to optimality; however, it fails to even find any feasible solution for . As before, these results are mostly explained by the large number of SDDP subproblems, and for these large-size instances,  and  have on average $114$ subproblems,  has $249$, and  has $593$.

Performance of the SDDP integrated B&C algorithm in large-size instances ($5\times6$ grid size).

4c|# Optimal     4c|# Feasible     4cOpt. Gaps (%)

Modality     1l|Cap.     1c     1c     1c     1c|     1c     1c     1c     1c|     1c     1c     1c     1c

3[2]*Type-1     1l|20%     2         0         0         0         8         10        7         1         53.4      51.8      56.1      82.0

1l|25%     5         2         0         0         5         8         6         0         35.7      26.2      27.9      -

1l|30%     10        8         0         0         0         1         4         0         -         15.4      4.3       -

3[2]*Type-2     1l|20%     3         0         0         0         7         10        7         1         43.5      52.9      61.4      86.0

1l|25%     6         3         0         0         4         7         6         0         37.2      24.9      32.8      -

1l|30%     9         8         0         0         1         1         6         0         12.0      13.3      6.4       -

2c|Total/Av.     35        21        0         0         25        37        36        2         36.4      30.7      31.5      84.0

§.§ Results for Approximation Methods

Next, we compare approximation methods to obtain feasible solutions for our aggregated models using 2SLDR and SDDP-based algorithms. We focus on the  transformation  because, as shown in Section <ref>, it achieves the best trade-off between policy quality and computational effort. Results for the other transformations can be found in Appendix <ref>.  We consider the 2SLDR variants described in Section <ref> (i.e., , , and )
and two SDDP-based methods to obtain a lower bound () and an upper bound ().  corresponds to the SDDP lower bounding procedure described in Section <ref>, and  takes the best integer solution found by  and performs additional calls to the SDDP sub-routine to obtain its exact objective value.

We first present results obtained by these approximation methods over mid-size instances in Table <ref>. We first analyze the average time to solve 2SLDR using Benders decomposition, and the average time to obtain a lower bound with . We see that  and  are the fastest among all options, and both take less time than . We also note that  takes significantly longer computational time than the other LDR variants, which is mostly explained by the large number of variables and distinct cost-to-go functions (see Section <ref> for a detailed explanation). Lastly,   is much faster than its exact variant , but it is still quite expensive compared to  and . Also, we note that  (i.e., evaluating the best integer solution of ) can take a significant computational effort due to the additional calls to the SDDP subroutine: around one to two hours for mid-size instances and up to $12$ hours for large-size instances.

Solution time and  quality of 2SLDR and SDDP bounds. Results for $\pa$ over mid-size instances.

4c|Average Time (sec)     5cRelative Difference (%)

Modality     1l|Cap.     1c     1c     1c|     1c|     1c     1c     1c     1c|     1c

3[2]*Type-1     1l|20%     813.0     164.4     157.9     549.1     0.12      0.26      0.11      0.01     0.36

1l|25%     380.1     84.3      77.7     347.0     0.08      0.11      0.03     0.26      0.71

1l|30%     408.2     79.0     85.9      400.7     0.12      0.25      0.00     0.55      1.24

3[2]*Type-2     1l|20%     1406.0     247.2     265.2     940.5     0.13      0.26      0.13      0.00     0.28

1l|25%     670.5     108.2     98.1     534.4     0.09      0.13      0.04     0.23      0.49

1l|30%     532.6     90.0      89.5     434.8     0.12      0.26      0.01     0.98      1.29

2c|Average     701.7     128.9     129.1     534.4     0.11      0.21      0.05     0.34      0.73

Table <ref> also shows the relative difference in the objective function for each approach with respect to the optimal solution found by , i.e., $(|Obj_{\extensive} - Obj_{i}|)/Obj_{\extensive}$ for all $i \in \{\thldr,\tldr,\mldr,\sddpub,\sddplb\}$. All the approximation methods obtain feasible solutions that are close to the true optimal solutions (optimality gaps are within 1%) and  obtains the best solutions in most instance configurations, except when the initial capacity level is low. Also,  and  dominate  in terms of solution quality and somewhat surprisingly,  obtains better solutions than  in most cases. Finally,  provides very strong lower bounds, making it a reliable lower-bounding approach when used together with LDR-based methods.

Solution time and quality of 2SLDR and SDDP bounds. Results for $\pa$ over large-size instances.

4c|Average Time (sec)     3cOpt. Gap (%)

Modality     Cap.      1c     1c     1c     (opt)     1c     1c     1c

3[2]*Type-1     1l|20%        5,353        3,732       13,564      (9)       0.62      0.39      0.24

1l|25%        4,235        4,716        13,557      (6)       2.56      2.06      1.91

1l|30%        1,555        1,675          2,067      (10)      0.80      0.63      0.27

3[2]*Type-2     1l|20%        7,489         4,855       18,323      (5)       9.00      7.33     9.01

1l|25%        6,664         5,560       14,790      (4)       4.98      4.38      4.28

1l|30%        1,717         1,640         6,971      (10)      0.81      0.63      0.30

2c|Av. (Total)        4,502         3,696       11,545      (44)      3.13      2.57     2.67

We now present results obtained by these approximation methods over large-size instances in Table <ref>. The LDR variant  is not included here as the corresponding formulation cannot be loaded into CPLEX due to the large number of first-stage LDR variables and constraints. Also, due to the lack of exact optimal solutions for large-size instances, the optimality gaps are computed with respect to the best lower bound found by , that is, we use $Obj_\sddplb$ instead of $Obj_\extensive$.

As in the mid-size instances, the 2SLDR variants can find high-quality solutions with relatively low computational effort. In fact,  and  find optimal solutions by the time limit for all the instances,  while  is able to prove optimality in 44 of the 60 tested instances (as indicated in the “(opt)" column). Nonetheless,  finds good lower bounds for most instances, and  yields better solutions in most instance configurations. These results suggest that 2SLDR, specifically , finds the best solutions with limited computational effort, while the SDDP-based methods can compute high-quality lower and upper bounds given enough computational effort.

§ CONCLUSIONS

We present a novel aggregation framework to address MSILP problems with mixed-integer state variables and continuous local variables. The main idea of the framework is to reformulate the original MSILP such that all integer state variables are first-stage and then apply a suitable transformation to aggregate these variables to reduce the problem size while not sacrificing the solution quality. We describe a number of transformations that rely on the underlying structure of the stochastic process, which we assume is given by an MC. Other transformations are possible in other problem settings (e.g., stage-independent or MC approximations of stage-dependent processes).

We propose a novel B&C framework integrated with the SDDP algorithm to solve the resulting aggregated MSILP. Also, we propose a computationally more tractable alternative based on the 2SLDR approximations and introduce new MC-based decision rules. These approaches are tested on a hurricane disaster relief logistics application with contingency modality decisions. Empirical results show the trade-off between different transformations regarding the solution quality and computational effort, and the best trade-off is yielded by transformations that consider information from current and previous MC states (i.e., ). The SDDP lower and upper bounding procedures have the best performance in terms of solution quality for large-scale instances, while the MC-based 2SLDR approach returns high-quality solutions with significantly less computational effort.

As shown by our empirical results, the performance of our aggregation framework in terms of solution quality and computational effort largely depends on the chosen transformation. Thus, a decision maker might need to try different alternatives to find a suitable transformation for their specific application. With this in mind, an interesting future research direction is to investigate how to choose an appropriate transformation for a given problem instance, possibly with the help of certain machine learning techniques. An alternative research direction is to develop an optimization framework to find the aggregation that has the fewest number of aggregated variables among all that lead to optimal policies.

This work was supported by Natural Sciences and Engineering Research Council of Canada [Grant RGPIN-2018-04984], the National Science Foundation of USA [Grant CMMI-2045744], and the National Center for Artificial Intelligence CENIA FB210017, Basal ANID Chile. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. Computations were performed on the Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto.

[Ahmed et al., 2020]
Ahmed S, Cabral FG, Freitas Paulo da Costa B (2020) Stochastic lipschitz
dynamic programming. Mathematical Programming 1–39.

[Alem et al., 2016]
Alem D, Clark A, Moreno A (2016) Stochastic network models for logistics
planning in disaster relief. European Journal of Operational Research

[Bally  Pagès, 2003]
Bally V, Pagès G (2003) A quantization algorithm for solving
multidimensional discrete-time optimal stopping problems. Bernoulli

[Bampou  Kuhn, 2011]
Bampou D, Kuhn D (2011) Scenario-free stochastic programming with polynomial
decision rules. 2011 50th IEEE Conference on Decision and Control and
European Control Conference, 7806–7812 (IEEE).

[Basciftci et al., 2019]
Basciftci B, Ahmed S, Gebraeel N (2019) Adaptive two-stage stochastic
programming with an application to capacity expansion planning. arXiv
preprint arXiv:1906.03513 .

[Ben-Tal et al., 2004]
Ben-Tal A, Goryashko A, Guslitzer E, Nemirovski A (2004) Adjustable robust
solutions of uncertain linear programs. Mathematical programming

[Birge, 1985]
Birge JR (1985) Decomposition and partitioning methods for multistage
stochastic linear programs. Operations research 33(5):989–1007.

[Birge  Louveaux, 1988]
Birge JR, Louveaux FV (1988) A multicut algorithm for two-stage stochastic
linear programs. European Journal of Operational Research

[Bodur et al., 2017]
Bodur M, Dash S, Günlük O, Luedtke J (2017) Strengthened benders cuts
for stochastic integer programs with continuous recourse. INFORMS
Journal on Computing 29(1):77–91.

[Bodur  Luedtke, 2022]
Bodur M, Luedtke JR (2022) Two-stage linear decision rules for multi-stage
stochastic programming. Mathematical Programming 191:347––380.

[Bonnans et al., 2012]
Bonnans JF, Cen Z, Christel T (2012) Energy contracts management by stochastic
programming techniques. Annals of Operations Research 200(1):199–222.

[Chen et al., 2008]
Chen X, Sim M, Sun P, Zhang J (2008) A linear decision-based approximation
approach to stochastic programming. Operations Research

[Davis et al., 2013]
Davis LB, Samanlioglu F, Qu X, Root S (2013) Inventory planning and
coordination in disaster relief efforts. International Journal of
Production Economics 141(2):561–573.

[Dowson, 2020]
Dowson O (2020) The policy graph decomposition of multistage stochastic
programming problems. Networks 76(1):3–23.

[Duran et al., 2011]
Duran S, Gutierrez MA, Keskinocak P (2011) Pre-positioning of emergency items
for care international. Interfaces 41(3):223–237.

[Füllner  Rebennack, 2021]
Füllner C, Rebennack S (2021) Stochastic dual dynamic programming and its
variants, preprint, Karlsruhe Institute of Technology.

[Füllner  Rebennack, 2022]
Füllner C, Rebennack S (2022) Non-convex nested benders decomposition.
Mathematical Programming 1–38.

[Garstka  Wets, 1974]
Garstka SJ, Wets RJB (1974) On decision rules in stochastic programming.
Mathematical Programming 7(1):117–143.

[Graumann et al., 2006]
Graumann A, Houston TG, Lawrimore JH, Levinson DH, Lott N, McCown S, Stephens
S, Wuertz DB (2006) Hurricane katrina: A climatological perspective:
Preliminary report .

[Hjelmeland et al., 2018]
Hjelmeland MN, Zou J, Helseth A, Ahmed S (2018) Nonconvex medium-term
hydropower scheduling by stochastic dual dynamic integer programming.
IEEE Transactions on Sustainable Energy 10(1):481–490.

[Kuhn et al., 2011]
Kuhn D, Wiesemann W, Georghiou A (2011) Primal and dual linear decision rules
in stochastic and robust optimization. Mathematical Programming

[Lodree Jr et al., 2012]
Lodree Jr EJ, Ballard KN, Song CH (2012) Pre-positioning hurricane supplies in
a commercial supply chain. Socio-Economic Planning Sciences

[Löhndorf  Shapiro, 2019]
Löhndorf N, Shapiro A (2019) Modeling time-dependent randomness in
stochastic dual dynamic programming. European Journal of Operational
Research 273(2):650–661.

[Morton, 1996]
Morton DP (1996) An enhanced decomposition algorithm for multistage stochastic
hydroelectric scheduling. Annals of Operations Research

[Pacheco  Batta, 2016]
Pacheco GG, Batta R (2016) Forecast-driven model for prepositioning supplies in
preparation for a foreseen hurricane. Journal of the Operational
Research Society 67(1):98–113.

[Pagès et al., 2004]
Pagès G, Pham H, Printems J (2004) An optimal markovian quantization
algorithm for multi-dimensional stochastic control problems.
Stochastics and dynamics 4(04):501–545.

[Pereira  Pinto, 1991]
Pereira MV, Pinto LM (1991) Multi-stage stochastic optimization applied to
energy planning. Mathematical programming 52(1):359–375.

[Philpott  De Matos, 2012]
Philpott AB, De Matos VL (2012) Dynamic sampling algorithms for multi-stage
stochastic programs with risk aversion. European Journal of operational
research 218(2):470–483.

[Philpott  Guan, 2008]
Philpott AB, Guan Z (2008) On the convergence of stochastic dual dynamic
programming and related methods. Operations Research Letters

[Regnier  MacKenzie, 2019]
Regnier ED, MacKenzie CA (2019) The hurricane decision simulator: A tool for
marine forces in new orleans to practice operations management in advance of
a hurricane. Manufacturing & Service Operations Management

[Sabbaghtorkan et al., 2020]
Sabbaghtorkan M, Batta R, He Q (2020) Prepositioning of assets and supplies in
disaster operations management: Review and research gap identification.
European Journal of Operational Research 284(1):1–19.

[Seraphin, 2019]
Seraphin H (2019) Natural disaster and destination management: The case of the
caribbean and hurricane irma. Current Issues in Tourism 22(1):21–28.

[Shapiro et al., 2021]
Shapiro A, Dentcheva D, Ruszczynski A (2021) Lectures on stochastic
programming: modeling and theory (SIAM).

[Shapiro  Nemirovski, 2005]
Shapiro A, Nemirovski A (2005) On complexity of stochastic programming
problems. Continuous optimization, 111–146 (Springer).

[Siddig  Song, 2022]
Siddig M, Song Y (2022) Multi-stage stochastic programming methods for adaptive
disaster relief logistics planning. arXiv preprint arXiv:2201.10678 .

[Taskin  Lodree Jr, 2010]
Taskin S, Lodree Jr EJ (2010) Inventory decisions for emergency supplies based
on hurricane count predictions. International Journal of Production
Economics 126(1):66–75.

[Van Slyke  Wets, 1969]
Van Slyke RM, Wets R (1969) L-shaped linear programs with applications to
optimal control and stochastic programming. SIAM journal on applied
mathematics 17(4):638–663.

[Yang et al., 2022]
Yang H, Duque D, Morton DP (2022) Optimizing diesel fuel supply chain
operations to mitigate power outages for hurricane relief. IISE
Transactions 1–14.

[Zhang  Sun, 2019]
Zhang S, Sun XA (2019) Stochastic dual dynamic programming for multistage
stochastic mixed-integer nonlinear optimization. arXiv preprint
arXiv:1912.13278 .

[Zou et al., 2018]
Zou J, Ahmed S, Sun XA (2018) Multistage stochastic unit commitment using
stochastic dual dynamic integer programming. IEEE transactions on Power
Systems 34(3):1814–1823.

[Zou et al., 2019]
Zou J, Ahmed S, Sun XA (2019) Stochastic dual dynamic integer programming.
Mathematical Programming 175(1):461–502.

[Zverovich et al., 2012]
Zverovich V, Fábián CI, Ellison EF, Mitra G (2012) A computational
study of a solver system for processing two-stage stochastic LPs with
enhanced benders decomposition. Mathematical Programming Computation

§ BENDERS DECOMPOSITION FOR THE  2SLDR MODEL

In this section we present a Benders decomposition approach to solve (<ref>). We consider a Benders master problem that handles the first-stage decisions and decompose the second-stage problem into one Benders subproblem per node. Because the number of subproblems can be exponentially many (e.g., when  is used), we consider an aggregated approach where we have one variable $\theta_{t,m}^L$ that approximates all the cost-to-go functions associated with nodes in stage $t$ and with MC state $m \in \markov_t$. This variant can be seen as a hybrid between the single-cut version (i.e., using a single $\theta$ variable) and the multi-cut version (i.e., using one $\theta_{t,n}^L$ variable for each subproblem) of Benders decomposition [Van Slyke  Wets, 1969, Birge  Louveaux, 1988]. Our preliminary experiments show that our Benders variant performs the best on our test instances, however, other variants of Benders decomposition and computational enhancements can also be considered (see, e.g., [Zverovich et al., 2012] and [Bodur et al., 2017]).

In what follows we present the details of our Benders decomposition for  (the other two cases,  and   are similar). The Benders master problem is given by

\begin{align}
\Qldrsub_{\rootnode} = \min\; &
\sum_{n\in \nodes}\prob_n \left( c_n^\top \zphi + d_n^\top(\varldr_{t(n)}^\top\xi^{t(n)}_n)  \right) + h_\rootnode^\top \yroot  + \sum_{t = 2}^\periods\sum_{m \in \markov_t}\theta_{t,m}^L \nonumber\\
\text{s.t.}\; &  \eqref{eq:ldr-constr0}-\eqref{eq:ldr-constr2}, \nonumber \\
& \theta^L_{t,m} \geq \alpha^L_{t,m}{}^\top \mu + \beta^L_{t,m}{}^\top \za + \gamma^L_{t,m}, \quad \forall  (\alpha^L_{t,m},\beta^L_{t,m},\gamma^L_{t,m})\in \coeffsldr_{t,m}, \ m \in  \markov_t, \ t \in \{2,...,\periods\}, \label{eq:ldr_benderscut}\\
& (\yroot, \varldr, \za)  \in  \R^r \times \R^{k\times \sum_{t\in \horizon}\dimxitt} \times \Z^{\ell \cdot \sum_{t \in \horizon} q_t}, \nonumber
\end{align}

where $\theta^L \in \R^{\sum_{t=2}^{\periods}|\markov_t|}$ is introduced to represent the cost-to-go function approximation for each Markovian state in all remaining stages (except for the first stage) via constraints (<ref>). Here $\coeffsldr_{t,m}$ stores the coefficients associated with the Benders cuts for each MC state $m \in \markov_t$.

1.2Benders Decomposition
[1]
Initialize B&B search for $\Qldrsub_\rootnode$ (master problem)
Choose a node B&B from its node list. Solve the associated node relaxation problem.
Current B&B node is infeasible or the relaxation bound is worse than the bound given by an incumbent solution
Prune this B&B node.
Current solution $(\varldrsol, \zasol, \thetasol^L)$ is integer
%% Cutting plane subroutine
$m \in \markov_t$ and $t \in \{2,...,\periods\}$
$n \in \staget(m)$
Solve $\Qldrsol_n:=\Qldr_n(\varldrsol, \zasol)$ and save dual solutions ${\pi^{n}}$.
Save value $\Qldrsol_{t,m}:=\sum_{n \in \staget(m)} p_n\Qldrsol_n$.
$|\Qldrsol_{t,m} - \thetasol^L_{t,m}| \geq \epsilon |\Qldrsol_{t,m}|$  %% Add cut if needed
Add inequality (<ref>) and exit the for loop.
No new inequality is added to the master problem $\Qldrsub_\rootnode$.
Branch on the current  B&B node and update list of nodes.
No more nodes in the B&B node list.

Algorithm <ref> describes our Benders decomposition for (<ref>). The procedure is a B&C algorithm similar to the one described in Algorithm <ref>. Once an integer solution is found during B&B search, we enter the cutting plane subroutine. The algorithm iterates over each stage and MC state and its corresponding nodes, where $\staget(m)$ is the set of all nodes in stage $t$ associated with MC state $m \in \markov_t$. It then evaluates the Benders subproblems for all the nodes associated with a particular $m$ and $t$ and checks if the current approximation $\thetasol^L_{t,m}$ is close enough to the actual value $\Qldrsol_{t,m}$ for a given tolerance $\epsilon$. If the approximation is not good enough, then the algorithm returns a cut and stops evaluating other subproblems. Similar to our SDDP variant, we add cuts to each candidate integer solution until no more cuts are found. The complete procedure ends when there are no more nodes to explore in the B&B search tree, i.e., we have found an optimal solution or proven infeasibility.

The cuts added to the master problem (<ref>) are Benders cuts over all the nodes associated to a given stage and MC state. Specifically, the cut coefficients for $t \in \{2,...,T\}$ and $m \in \markov_t$ are given by

\begin{align*}
\alpha_m^L & = \sum_{n \in \staget(m)}\prob_n {\pi^{n}}^\top  A_{n}, \qquad \beta_m =  \sum_{n \in \staget(m)} \prob_n {\pi^{n}}^\top B_{n}, \\
\gamma_m  &= \sum_{n \in \staget(m)} \prob_n \left( \Qldrsol_n - ({\pi^{n}}^\top  A_{n} )^\top \hat{\mu}_{t(n) -1}- ({\pi^{n}}^\top B_{n})^\top \zasol_{\phi_t(a(n))} \right),
\end{align*}

where ${\pi^{n}}$ are the dual solutions associated with the cost-to-go function $\Qldr(n)$. Note that each coefficient here corresponds to a weighted sum of the associated coefficients in the subproblems. Also, the algorithm generates Benders feasibility cuts when necessary, which is common in 2SLDR approximations [Bodur  Luedtke, 2022].

§ PROBLEM DESCRIPTION AND INSTANCE GENERATION FOR HDR

We now describe the main components of our HDR problem and present details of the instance generation procedure. We  focus on the grid representation for the problem, the MC, and three main parameters (i.e., demand, modalities, and capacity). All other details about instances generation (e.g., cost parameters) can be found in our instances generation code. We will make the instance generation code and the set of instances publicly available upon publication.

§.§ Grid Representation and Locations of Shelters and DCs

Similar to previous works in the literature, we use a grid network to represent the potential locations of shelters and DCs, as well as possible locations of the hurricane (see Figure <ref>). The top row on the grid represents the land, where each cell contains a subset of shelters and DCs. The remaining cells in the grid correspond to the sea and are used to represent possible locations of the hurricane at any stage during the planning horizon. For simplicity, we assume that any shelter can be supplied by any DC in the network.

A $6\times 4$ grid example that shows the initial location of the hurricane $(m^x,m^y)=(5,0)$ and all possible locations in the following stages (i.e., shaded region).

Land and sea cells have pre-defined dimensions: land cells have a size of $50\times100$, while sea cells have a size of $20\times100$. We randomly generate 3 to 7 shelters and 2 to 4 DCs in each land cell. Shelters and DCs are positioned uniformly at random inside their corresponding cells. The coordinates of the hurricane located in a cell are assumed to be the coordinates of the center of the cell.

§.§ MC description

As described in the main text, each MC state corresponds to two attributes of the hurricane: location and intensity, that is, $m=(m^x, m^y, m^i)\in \markov$ where $m^x$ and $m^y$ are the $x$- and $y$-coordinates of the hurricane's location, respectively, and $m^i$ represents the hurricane's intensity level. The MC is given by two independent transition probability matrices: one for the hurricane intensity level and one for its location. The intensity level transition matrix is identical to the one used by [Pacheco  Batta, 2016], which considers six levels of intensity $\{0,1,\dots,5\}$ where level 5 corresponds to the maximum hurricane intensity and level 0 corresponds to the case when the hurricane dissipates. The following matrix shows the transition probabilities where $P_{ij}$ represents the probability of transitioning from intensity level $i$ to $j$.

\[
P = \bordermatrix{
~ & 0   & 1   & 2   & 3   & 4   & 5     \cr
0 & 1   & 0   & 0   & 0   & 0   & 0     \cr
1 & 0.11& 0.83& 0.06& 0   & 0.6 & 0     \cr
2 & 0   & 0.15& 0.6 & 0.25& 0   & 0     \cr
3 & 0   & 0   & 0.04& 0.68& 0.28& 0     \cr
4 & 0   & 1   & 0   & 0.18& 0.79& 0.03  \cr
5 & 0   & 0   & 0   & 0   & 0.5 & 0.5   \cr
\]