Generating Lagrangian Cuts Using Normalized Dual Problems in Multistage Stochastic Mixed-Integer Programming

1. Introduction

In this paper, we study the generation of Lagrangian cuts in decomposition methods for solving multistage stochastic mixed-integer linear programs (MS-MILPs), such as stochastic dual dynamic integer programming (SDDiP) (Zou et al. 2019). We present a framework to generate Lagrangian cuts with favorable properties by solving special normalized Lagrangian dual problems.

1.1. Motivation and Prior Work

Multistage stochastic programs are relevant to model decision-making processes in practice because often, sequential decisions have to be made over a finite number of stages and under uncertainty considering the problem data of the following stages. For the solution of these problems, tailored decomposition methods are most widely used; stochastic dual dynamic programming (SDDP) is among the most prominent ones (Pereira and Pinto 1991). These methods decompose the large-scale problem into stage- and scenario-specific subproblems coupled by state variables and so-called value functions, which are denoted $Q_n(·)$. For linear problems, these functions are convex polyhedral and can be exactly represented by finitely many affine functions called cuts (Benders 1962). However, in many applications, some of the decision variables have to be integer or binary. In this case, we obtain an MS-MILP, and the value functions are in general nonconvex and discontinuous.

A key challenge is that in this case, linear cuts, such as Benders cuts, are in general not sufficient to ensure (almost sure) finite convergence of SDDP-like decomposition methods as they are no tight approximators of $Q_n(·)$. Recently, more focus has been put on Lagrangian cuts, which are constructed by solving special Lagrangian dual problems and stronger than Benders cuts (Zou et al. 2019) as they can recover the closed convex envelope $\overline{\text{co}}(Q_n)(·)$ of $Q_n(·)$ (Füllner et al. 2024). Even more, in the special case that all state variables are binary, Lagrangian cuts are tight for $Q_n(·)$. This is at the core of the SDDiP algorithm. Using a binary approximation of the state space, SDDiP can also be applied to general MS-MILPs with convergence guarantees (Zou et al. 2019).

Nonetheless, applying Lagrangian cuts in practice comes with some considerable challenges.

1. Even if the Lagrangian dual is a convex optimization problem, it may be very costly to solve repeatedly because of its nonsmooth objective function. This is only aggravated if the state space is artificially increased by a binary approximation for general MS-MILPs.

2. Lagrangian dual problems are often degenerate with multiple optimal solutions, and the cuts associated with some of these solutions can be very weak. This issue is especially common for the binary state space required in SDDiP because all cuts are constructed at extreme points. It is, therefore, important not only to solve the dual problem but also, to identify a solution that yields a strong cut.

3. Whereas tight cuts are sufficient to ensure theoretical convergence of SDDiP, (solely) relying on them in practice may lead to very slow convergence. At worst, a complete enumeration of the binary state space is required. It is plausible that there exist alternative (even nontight) cuts that may significantly speed up the convergence process.

Some first attempts to address these challenges have been made recently, mostly with respect to (a) and (c). In the original SDDiP paper, it was proposed to combine tight cuts with strengthened Benders cuts that are not tight in general but stronger than Benders cuts and efficient to compute (Zou et al. 2019). Similarly, a scheme alternating between Lagrangian and Benders cuts is proposed by Bansal and Küçükyavuz (2026). Rahmaniani et al. (2020) present a heuristic to generate Lagrangian cuts more efficiently using inner approximations or partial relaxations. For two-stage stochastic problems, Chen and Luedtke (2022) suggest to restrict the dual feasible set to the span of Benders cut coefficients of previous iterations, thereby losing tightness guarantees but obtaining reasonably good Lagrangian cuts in short time. With respect to (b), Deng and Xie (2024) propose to first generate cheap and tight but weak integer optimality cuts and then, to solve some auxiliary linear programs (LPs) to strengthen them in order to obtain strong Lagrangian cuts.

In addition, there has been a lot of research on alternative cut-generation techniques addressing challenge (b), yet for Benders cuts in Benders decomposition (BD) instead of Lagrangian cuts, see Table 1. To address degeneracy and weak cuts, Magnanti and Wong (1981) present a two-step approach to generate Pareto-optimal cuts. Their ideas are improved by Papadakos (2008), who shows that it is sufficient to solve a single optimization problem. Sherali and Lunday (2013) propose to generate certain maximally nondominated Benders cuts by solving a perturbation of the original subproblem. Glomb et al. (2026) introduce so-called optimal line-shifting cuts satisfying a novel notion of Pareto optimality.

Table 1. Examination of Cut-Quality Criteria in the Literature on BD and Disjunctive Programming

A novel epigraph-based perspective on generating Benders cuts is introduced in Fischetti et al. (2010). In contrast to standard BD, by using a formulation with an additional dual multiplier, the proposed cut-generation problem allows us to generate optimality and feasibility cuts in a unified way. On the flip side, the cut-generation problem is unbounded without further adjustments. Therefore, the authors suggest to solve a special normalized version instead.

The formulation from Fischetti et al. (2010) allows for different normalization techniques and by that, the generation of cuts satisfying different quality criteria. Hosseini and Turner (2025) show how normalizations can be used to generate deep cuts in BD, which are characterized by maximizing the distance between the separating hyperplane and the incumbent to separate, a property that is also explored in Cornuéjols and Lemaréchal (2006) and Cadoux (2010). By relating the cut generation to the so-called reverse polar set, Brandenberg and Stursberg (2021) show how facet-defining and Pareto-optimal cuts can be generated in BD using linear normalizations (LNs). The same cut-generation procedure is put forward by Seo et al. (2022) but motivated from a different geometric perspective. An alternative yet related approach to generate facet-defining cuts is presented by Conforti and Wolsey (2019).

Many of these approaches are reported to lead to significant performance improvements compared with standard BD. Therefore, it seems reasonable to explore how they (or similar ideas) can be translated to stochastic programming and Lagrangian cuts. A first attempt in this direction is the work by Chen and Luedtke (2022). They apply the unified cut-generation perspective from Fischetti et al. (2010) and a special normalization to generate Lagrangian cuts in branch-and-cut methods to solve two-stage stochastic mixed-integer linear programs (MILPs).

1.2. Contribution

In this paper, similar to Chen and Luedtke (2022), we apply the unified perspective from Fischetti et al. (2010) to the multistage stochastic case, and based on this, we present a general framework with which Lagrangian cuts of various styles can be generated within decomposition methods for MS-MILPs. We show that as for BD, different normalization techniques in the Lagrangian dual lead to Lagrangian cuts with distinct properties, such as maximal depth, Pareto optimality, or being facet defining with respect to $\overline{\text{co}}(Q_n)(·)$.

Compared with the work by Chen and Luedtke (2022), our work differs in several aspects. First, we consider a formulation for the multistage case, even though the extension from the two-stage case is straightforward. Second, we consider a variety of possible normalization techniques for the Lagrangian dual, including different norms and linear functions. In contrast, Chen and Luedtke (2022) focus on a specific type of normalization using the 1-norm that is directly linked to their idea of restricting the dual feasible set. Third, they solve two-stage stochastic problems using a branch-and-cut method. Stronger cuts can help accelerate the solution process in these methods but are not required to achieve convergence. Therefore, the authors’ main focus is the computationally efficient generation of reasonably strong (yet not necessarily supporting) Lagrangian cuts, thus addressing challenges (a) and (c). We, on the other hand, explore in detail the effect of different normalizations on the theoretical properties of the Lagrangian cuts with the aim of satisfying certain quality criteria. Hence, our work mostly addresses challenge (b) above. As the new Lagrangian cuts do not satisfy the traditional notion of tightness but lead to better lower bounds (LBs) in our computational tests, our work also refers to challenge (c). Even more, we show that using our proposed cut-generation framework, convergence of SDDiP can still be assured.

Compared with the existing literature on BD, our work differs as well. First, we apply the cut-generation formulation from Fischetti et al. (2010) to the stochastic and Lagrangian setting. Second, we apply various normalization concepts from the literature on BD and show how they and the properties of the obtained cuts translate to this setting. A notable difference in the Lagrangian setting is that we do not tightly approximate $\text{epi}(Q_n)$ itself but $\text{epi}(\overline{\text{co}}(Q_n))$ for which no explicit description is readily available. This also limits the applicability of certain techniques that are tailored to BD. Finally, the MS-MILP setting poses several challenges for the computation of core points that the linear normalizations that we use rely on, among them infeasibility as well as aggravated verification of core-point candidates. For this reason, we propose several heuristics to identify core points in our specific setting. We run extensive computational tests to evaluate their performance.

The key contributions of this paper are summarized again below.

1. We show how the unified cut-generation formulation for BD from Fischetti et al. (2010) can be applied to the generation of Lagrangian cuts for MS-MILPs. This idea has already been used by Chen and Luedtke (2022) in two-stage stochastic MILPs, but it has not been extended to the multistage case yet.

2. As the obtained Lagrangian dual problems are unbounded, some normalization is required to select cut coefficients in a reasonable way. We draw on recent concepts for cut selection in BD, such as optimizing over the reverse polar set (Brandenberg and Stursberg 2021) or generating deep cuts (Hosseini and Turner 2025), and we translate them to the stochastic and Lagrangian setting. This way, we obtain a variety of different normalization techniques and by that, generalize the cut-generation approach from Chen and Luedtke (2022). We further show that different normalization techniques in the Lagrangian dual lead to Lagrangian cuts with distinct properties, such as maximal depth, Pareto optimality, or being facet defining with respect to $\overline{\text{co}}(Q_n)(·)$.

3. We show that linear normalizations are closely related to the identification of core points in the epigraphs $\text{epi}(Q_n)$, which can be challenging for MS-MILPs. Therefore, we propose four heuristic approaches for the computation of core-point candidates and study them computationally.

4. Our proposed cut-generation framework can be incorporated into SDDiP and related decomposition methods. We prove that under some assumptions, still (almost sure) finite convergence of these methods can be guaranteed.

5. We perform extensive computational tests for SDDiP incorporating our cut-generation approach on two known test problems from the literature. We show that the obtained lower bounds in SDDiP are majorly improved compared with conventional cuts. However, we also observe that this does not in all cases lead to an improvement of the obtained policies and that solving the Lagrangian dual problem remains a computational challenge.

For reasons of space, some details, including all technical proofs, are moved to the Online Appendix.

2. Problem Formulation

We start by introducing MS-MILPs and their decomposition formally. We consider a finite number $T\in \mathbb{N}$ of stages, where some of the problem data are uncertain and evolve according to a known stochastic process $\xi ≔({\xi }_1,\dots ,{\xi }_T)$ with deterministic ${\xi }_1$. We assume that the random data vectors ${\xi }_t,t=1,\dots ,T,$ are discrete and finite such that the uncertainty can be modeled by a finite scenario tree. Let $\mathcal{N}$ denote the set of nodes of this tree. For each node $n\in \mathcal{N}$, the unique ancestor node is denoted by $a(n)$, and the set of child nodes is denoted by $\mathcal{C}(n)$. The probability for some node n is $p_n>0$ and assumed to be known. The transition probabilities between adjacent nodes $n,m\in \mathcal{N}$ can then be determined as $p_{nm}≔\frac{p_m}{p_n}$. For the root node r, we assume $a(r)=\varnothing$ and $p_r=1$. We define $\overline{\mathcal{N}}≔\mathcal{N}\setminus \{r\}$ to address the set of nodes without the root node and $\tilde{\mathcal{N}}$ to address the set of nodes without leaf nodes, and we denote by $\mathcal{N}(t)$ the nodes at stage t.

For each node $n\in \mathcal{N}$, we distinguish state variables $x_n\in {\mathbb{R}}^{d_n}$, which also appear in child nodes of n, and local variables $y_n\in {\mathbb{R}}^{{\tilde{d}}_n}$. $f_n(·)$ denotes the objective function of node n, and ${\mathcal{F}}_n(x_{a(n)})$ denotes the feasible set of node n, which depends on the state variable $x_{a(n)}$ from the ancestor node. We assume that $f_n(·)$ is a linear function in $x_n$ and $y_n$ and that ${\mathcal{F}}_n(·)$ is a mixed-integer polyhedral set for all $x_{a(n)}$ defined by

Here, $A_n,B_n,C_n,b_n$ denote appropriately defined data matrices and vectors. The sets $X_n$ and $Y_n$ are intersections of polyhedral sets ${\overline{X}}_n,{\overline{Y}}_n$ and possible integrality constraints. In the following, we also refer to $X_n$ as the state space.

An MS-MILP can then be expressed by its dynamic programming equations. For the root node, we obtain

with $x_{a(r)}=0$, and $v^{*}$ is the optimal value of the original problem. Let $\overline{\mathbb{R}}≔\mathbb{R}\cup \{+\infty \}$. For all $n\in \overline{\mathcal{N}}$, the value function $Q_n:{\mathbb{R}}^{d_{a(n)}}\to \overline{\mathbb{R}}$ is defined by

For the leaf nodes $n\in \mathcal{N}\setminus \tilde{\mathcal{N}}$, we set ${\sum }_{m\in \mathcal{C}(n)}{p}_{nm}{Q}_m(x_n)\equiv 0$. Moreover, we set $Q_n(x_{a(n)})=+\infty$ if ${\mathcal{F}}_n(x_{a(n)})=\varnothing$, and we denote by $\text{dom}(Q_n)$ the effective domain of $Q_n(·)$. Furthermore, we denote by $\overline{\text{co}}(Q_n)(·)$ its closed convex envelope, the point-wise supremum of all affine functions majorized by $Q_n(·)$.

Remark 1.

Note that regarding $Q_n(·)$ as a function on ${\mathbb{R}}^{d_{a(n)}}$ is not standard in stochastic programming. Often, it is (implicitly) assumed to be defined only on the domain $X_{a(n)}$. However, from our view, allowing $Q_n(·)$ to be defined on ${\mathbb{R}}^{d_{a(n)}}$ with extended real values appears more suitable for the following steps.

As this proves beneficial in the cut-generation process, we introduce local variables $z_n$ and accompany them with copy constraints $x_{a(n)}=z_n$ and constraints $z_n\in Z_{a(n)}$, with $Z_{a(n)}\supseteq X_{a(n)}$. The most natural choice is $Z_{a(n)}=X_{a(n)}$, but also, other choices are possible (e.g., $Z_{a(n)}=\text{conv}(X_{a(n)})$). For more details, we refer to Füllner et al. (2024). This reformulation yields the equivalent subproblems

where ${\mathcal{F}}_n≔\{(x_n,y_n,z_n)\in {\mathbb{R}}^{d_{a(n)}}\times {\mathbb{R}}^{d_n}\times {\mathbb{R}}^{{\tilde{d}}_n}\hspace{0.17em}:\hspace{0.17em}(x_n,y_n)\in {\mathcal{F}}_n(z_n)\}$.

For the remainder of this article, we make some basic assumptions.

Assumption 1.

The following conditions are satisfied by (1)–(3).

A1. For all $n\in \mathcal{N}$, the sets $X_n$ and $Y_n$ are compact.

A2. For all $n\in \mathcal{N}$, all coefficients in $A_n,B_n,C_n,b_n$, $f_n$, ${\overline{X}}_n$, and ${\overline{Y}}_n$ are rational.

A3. For all $n\in \overline{\mathcal{N}}$, $Z_{a(n)}$ is compact, rational MILP representable, and $Z_{a(n)}\subseteq \text{dom}(Q_n)$.

Note that property (A1) of Assumption 1 immediately implies that $F_n(x_{a(n)})$ is bounded for all $x_{a(n)}\in {\mathbb{R}}^{d_{a(n)}}$ and $n\in \mathcal{N}$. Property (A3) of Assumption 1 and the way that $Q_n(·)$ is defined imply relatively complete recourse, a standard assumption in stochastic programming. We obtain the following standard properties for the value functions.

Lemma 1.

For all $n\in \overline{\mathcal{N}}$, the value functions $Q_n(·)$ are proper, lower semicontinuous, and piecewise polyhedral with finitely many pieces.

By applying the properness reasoning to the root node, we conclude that $v^{*}$ is finite.

3. Generating Lagrangian Cuts with Normalized Dual Problems

In this section, we present a framework for the generation of Lagrangian cuts when solving MS-MILPs that serves as an alternative to the classical Lagrangian cut-generation framework from SDDiP (Zou et al. 2019) (see Online Appendix EC.1 for comparison). The framework is based on three key ingredients.

• It uses the idea of the unified cut-generation formulation for BD from Fischetti et al. (2010) and applies it to MS-MILPs to derive a special Lagrangian dual problem. This is the subject of Sections 3.1 and 3.2. The dual problem generalizes the one used by Chen and Luedtke (2022) to the multistage case.

• As this dual problem is unbounded, our framework also follows the aforementioned works in applying a normalization. Deviating from the static normalization ${\sum }_{k=0}^{d_n}{\pi }_{nk}=1$ (Fischetti et al. 2010) and the ${\mathcal{\ell }}^1$-norm normalization (Chen and Luedtke 2022) used in these works, our framework allows for general normalizations. We discuss different options using norms and linear functions, and we explore the impact on the properties of the obtained Lagrangian cuts. This is the subject of Sections 3.3–3.5.

• When a linear normalization is used, knowledge of a core point within $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ is required, which can be challenging to achieve for MS-MILPs. We thus present different heuristics for that purpose. This is the subject of Section 3.6.

Algorithm 1 provides a general overview on the cut-generation framework.

Algorithm 1

(Cut-Generation Framework for Lagrangian Cuts)

Require: Subproblem (6) for node $n\in \tilde{\mathcal{N}}$ with index i. NORM = true or NORM = false.

• 1: Solve Subproblem (6) to obtain $x_n^i$ and ${\theta }_m^i$ for all $m\in \mathcal{C}(n)$.

• 2: for child node $m\in \mathcal{C}(n)$ do

• 3:  ifNORM = true then

• 4:   Choose some norm $\Vert ·\Vert$, and set up the dual problem (11) for m and i with $g_m({\pi }_m,{\pi }_{m0})=\Vert {\pi }_m,{\pi }_{m0}\Vert$.

• 5:  else

• 6:   Use heuristics Mid, Eps, Relint, or Conv from Section 3.6 to compute a core-point candidate $({\check{x}}_n^i,{\check{\theta }}_m^i)$.

• 7:   Compute coefficients $(u_m,u_{m0})=({\check{x}}_n^i,{\check{\theta }}_m^i)-(x_n^i,{\theta }_m^i)$.

• 8:   Set up the dual problem (11) for m and i with $g_m({\pi }_m,{\pi }_{m0})=u_m^{\top }{\pi }_m+u_{m0}{\pi }_{m0}$.

• 9:  Solve the dual problem (11) for m and i, and store the multipliers $({\pi }_m^i,{\pi }_{m0}^i)$.

• 10:  Construct a cut of form (10).

• 11:  Use the cut to update the set ${\mathrm{\Psi }}_m^i$ to ${\mathrm{\Psi }}_m^{i+1}$ in Subproblem (6) for node n.

3.1. An Epigraph Perspective on Cut Generation

We first introduce the cut-generation perspective from Fischetti et al. (2010) yet tailor it to our case of MS-MILPs. Recall the definition of the epigraph of the value functions $Q_n(·),n\in \overline{\mathcal{N}}$:

We can use this definition to reformulate Subproblems (3) to

Here and in the following, we use $({\theta }_m)$ as a shortened notation for ${({\theta }_m)}_{m\in \mathcal{C}(n)}$.

Remark 2.

The condition $x_n\in \text{dom}(Q_m)$ from the definition in (4) is always satisfied implicitly in Problem (5) because $\text{dom}(Q_m)=Z_n$ by Assumption 1 and $x_n\in X_n\subseteq Z_n$ by construction.

In classical BD, iteratively optimality cuts are constructed to approximate $Q_n(·)$, and if required, feasibility cuts are constructed to approximate $\text{dom}(Q_n)$. This is done by solving distinct cut-generation problems. However, from (4), it is evident that both types of cuts actually approximate $\text{epi}(Q_n)$. Therefore, we may as well consider a unified cut-generation problem to obtain polyhedral approximations ${\mathrm{\Psi }}_m$ of the sets $\text{epi}(Q_m)$ (Fischetti et al. 2010). Although the need for feasibility cuts is ruled out by property (A3) of Assumption 1 in our setting here, the unified cut-generation framework still proves itself valuable as we shall see.

Remark 3.

In multistage stochastic programming, cuts are often aggregated to obtain single cuts for the expected value functions ${\sum }_{m\in \mathcal{C}(n)}{p}_{nm}{Q}_m(x_n)$ (single-cut approach) as this reduces the total number of cuts in the subproblems. The unified perspective here requires a separate set of approximations ${\mathrm{\Psi }}_m$ for each $\text{epi}(Q_m)$ (multicut approach) instead as it uses distinct values ${\theta }_m^i$ in the cut-generation process.

As the value functions $Q_n(·)$ are not known explicitly, for the cut-generation process, we replace each occurrence of $\text{epi}(Q_m)$ in (5) with its current approximation ${\mathrm{\Psi }}_m^{i+1}$ (with iteration index i). This way, we actually generate cuts for $\text{epi}({\underset{¯}{Q}}_m^{i+1})$. However, by construction, they also yield outer approximations of $\text{epi}(Q_m)$. To avoid unboundedness, each set is initialized with a valid outer approximation ${\mathrm{\Psi }}_m^0,m\in \mathcal{C}(n)$.

For notational simplicity, for the remainder of this paper, we define the set

and we further define ${\lambda }_n≔(x_n,y_n,({\theta }_m))$ and $c_n^{\top }{\lambda }_n≔f_n(x_n,y_n)+{\sum }_{m\in \mathcal{C}(n)}{p}_{nm}{\theta }_m$ (recall that $f(·)$ is linear). Then, the approximate subproblem associated with (5) can be compactly written as

We make another assumption for the remainder of this paper.

Assumption 2.

For all $n\in \overline{\mathcal{N}}$ and all iterations i, all linear cuts defining the polyhedral set ${\mathrm{\Psi }}_m^{i+1}$ are defined by rational coefficients.

Furthermore, we later require the following result for the convex hull $\text{conv}({\mathcal{W}}_n^{i+1})$ of set ${\mathcal{W}}_n^{i+1}$.

Remark 4.

The set $\text{conv}({\mathcal{W}}_n^{i+1})$ is a closed convex polyhedron. That means that there exist matrices ${\tilde{A}}_n,{\tilde{B}}_n$ and a vector ${\tilde{d}}_n$ such that $\text{conv}({\mathcal{W}}_n^{i+1})=\{({\lambda }_n,z_n)\hspace{0.17em}:\hspace{0.17em}{\tilde{A}}_n{\lambda }_n+{\tilde{B}}_n{z}_n\ge {\tilde{d}}_n\}$.

3.2. An Unbounded Lagrangian Dual Problem

We now address the cut-generation process in the unified framework. Given a point $(x_{a(n)}^i,{\theta }_n^i)$, we consider the feasibility problem

which can be shown to verify if $(x_{a(n)}^i,{\theta }_n^i)\in \text{epi}({\underset{¯}{Q}}_n^{i+1})$; see Online Appendix EC.2.1 for a formal proof.

To generate Lagrangian cuts, we apply a Lagrangian relaxation to Problem (7). A key difference to the classical Lagrangian relaxation from SDDiP (see Online Appendix EC.1) is that not only $x_{a(n)}^i$, but $(x_{a(n)}^i,{\theta }_n^i)$ is regarded as a fixed incumbent for the cut-generation process. Therefore, we relax all constraints containing either $x_{a(n)}^i$ or ${\theta }_n^i$. For given multipliers $({\pi }_n,{\pi }_{n0})\in {\mathbb{R}}^{d_{a(n)}}\times {\mathbb{R}}^{+}$, the dual function is then given by

The Lagrangian dual problem, which corresponds to the one used by Chen and Luedtke (2022) for the two-stage case, is

We state some important properties of this dual problem, with a proof provided in Online Appendix EC.2.2.

Theorem 1.

The optimal value ${\widehat{v}}_n^{D,i+1}(x_{a(n)}^i,{\theta }_n^i)$ of the Lagrangian dual (9) satisfies

$${\widehat{v}}_n^{D,i+1}(x_{a(n)}^i,{\theta }_n^i)=\{\begin{array}{ll}0,\hfill & \mathit{\text{if}} (x_{a(n)}^i,{\theta }_n^i)\in \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))\hfill \\ +\infty ,\hfill & \mathit{\text{else}}.\hfill \end{array}$$

Theorem 1 implies that the Lagrangian dual is unbounded whenever $(x_{a(n)}^i,{\theta }_n^i)\notin \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$. Therefore, there exists an unbounded ray $({\pi }_n^i,{\pi }_{n0}^i)$ such that ${ℒ}_n^{i+1}({\pi }_n^i,{\pi }_{n0}^i)-{({\pi }_n^i)}^{\top }{x}_{a(n)}^i-{\pi }_{n0}^i{\theta }_n^i>0$, and $(x_{a(n)}^i,{\theta }_n^i)$ violates the following Lagrangian cut.

Definition 1.

For all $n\in \overline{\mathcal{N}}$ and some multipliers $({\pi }_n^i,{\pi }_{n0}^i)$, a Lagrangian cut is given by

$${\pi }_{n0}^i{\theta }_n+{({\pi }_n^i)}^{\top }{x}_{a(n)}\ge {ℒ}_n^{i+1}({\pi }_n^i,{\pi }_{n0}^i).$$(10)

This cut is valid for any feasible $({\pi }_n^i,{\pi }_{n0}^i)$ in (9) as proven in Online Appendix EC.2.3.

Lemma 2.

For any $({\pi }_n^i,{\pi }_{n0}^i)\in {\mathbb{R}}^{d_{a(n)}}\times {\mathbb{R}}^{+}$, the Lagrangian cut (10) is satisfied by all $(x_{a(n)},{\theta }_n)\in \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ and thus, by all $(x_{a(n)},{\theta }_n)\in \text{epi}(Q_n)$.

3.3. A Normalized Lagrangian Dual Problem

It is not immediately clear how to select cut coefficients $({\pi }_n^i,{\pi }_{n0}^i)$ from the unbounded Lagrangian dual (9) in the most reasonable way. Therefore, we consider a normalized problem instead.

Definition 2.

For some homogeneous normalization function $g_n:{\mathbb{R}}^{d_n}\times {\mathbb{R}}^{+}\to \mathbb{R}$, the normalized Lagrangian dual is defined as

$$\begin{array}{l}{\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)≔\underset{{\pi }_n,{\pi }_{n0}}{\max }\{{ℒ}_n^{i+1}({\pi }_n,{\pi }_{n0})-{\pi }_n^{\top }{x}_{a(n)}^i-{\pi }_{n0}{\theta }_n^i\hspace{0.17em}:\hspace{0.17em}{g}_n({\pi }_n,{\pi }_{n0})\le 1,{\pi }_{n0}\ge 0\}.\hfill \end{array}$$(11)

Remark 5.

As long as the normalization constraint $g_n({\pi }_n,{\pi }_{n0})\le 1$ is satisfied by some neighborhood N of the origin, we do not exclude any potential cuts because of the scaling property of ${\pi }_{n0}$; see Remark EC.1 in Online Appendix EC.1 and Chen and Luedtke (2022).

Remark 6.

Fischetti et al. (2010) use $g_n({\pi }_n,{\pi }_{n0})={\sum }_{k=0}^{d_n}{\pi }_{nk}=1$ and an adjustment accounting for zero entries in the problem matrices. Chen and Luedtke (2022) propose $g_n({\pi }_n,{\pi }_{n0})={\Vert {\pi }_n,{\pi }_{n0}\Vert }_1$ with a possible adjustment tailored to their reduction of the dual space.

Remark 7.

Problem (11) can be solved using the Kelley cutting-plane method (Kelley 1960) as presented in Online Appendix EC.6 or using related bundle methods (Lemaréchal et al. 1995).

The normalized problem has the following important property.

Lemma 3.

If $(x_{a(n)}^i,{\theta }_n^i)\in \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$, then ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)=0$ and vice versa.

We provide a proof in Online Appendix EC.2.4. Lemma 3 allows us to solely focus on the case where $(x_{a(n)}^i,{\theta }_n^i)\notin \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ for the remainder of this section.

The reasons to consider a normalization of the unbounded Lagrangian dual (9) are twofold. On the one hand, we aspire to determine coefficients $({\pi }_n^i,{\pi }_{n0}^i)$ by solving a bounded and feasible optimization problem. A common approach to achieve this is to fix its unbounded objective to one. Combined with Remark 4, this allows us to identify unbounded rays by analyzing a compact polyhedron (Fischetti et al. 2010). Alternatively, it can be achieved by bounding Problem (9) artificially (e.g., by introducing bounds on the dual multipliers or by introducing a carefully designed normalizing constraint).

On the other hand, in the light of challenge (b) presented in Section 1.1, we want to make sure that the obtained cuts do not only separate the incumbent $(x_{a(n)}^i,{\theta }_n^i)$ from $\text{epi}(Q_n)$ but also, are of good approximation quality. As stated in Table 1 in Section 1.1, various quality criteria for cutting planes have been put forward in the literature. We focus on three important criteria here.

• Maximizing cut depth. Deep cuts realize a maximum distance between the incumbent $(x_{a(n)}^i,{\theta }_n^i)$ and the separating hyperplane (i.e., they cut as deeply as possible into the suboptimal region).

• Being facet defining. Facet-defining cuts reproduce (some of the finitely many) facets of a convex polyhedral set: in our case, $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$.

• Pareto optimality (Magnanti and Wong 1981). For ${\pi }_{n0}>0$, Pareto-optimal cuts (10) are nondominated in the sense that there exists no other cut ${\tilde{\pi }}_{n0}{\theta }_n+{({\tilde{\pi }}_n)}^{\top }{x}_{a(n)}\ge {ℒ}_n^{i+1}({\tilde{\pi }}_n,{\tilde{\pi }}_{n0})$ such that

  $$\frac{1}{{\tilde{\pi }}_{n0}}({ℒ}_n^{i+1}({\tilde{\pi }}_n,{\tilde{\pi }}_{n0})-{({\tilde{\pi }}_n)}^{\top }{x}_{a(n)})\ge \frac{1}{{\tilde{\pi }}_{n0}}({ℒ}_n^{i+1}({\pi }_n^i,{\pi }_{n0}^i)-{({\pi }_n^i)}^{\top }{x}_{a(n)})$$

for all $(x_{a(n)},{\theta }_n)\in \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$. Importantly, Pareto optimality with respect to $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ does not necessarily imply Pareto optimality with respect to $\text{epi}({\underset{¯}{Q}}_n^{i+1})$ but is easier to achieve.

As shown in the literature on BD, many of these criteria can be satisfied by optimizing over the so-called reverse polar set of $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ shifted by the incumbent $(x_{a(n)}^i,{\theta }_n^i)$, which provides a characterization of normal vectors of $(x_{a(n)}^i,{\theta }_n^i)$-separating hyperplanes and thus, is an important tool in the theory on cut generation (Cornuéjols and Lemaréchal 2006, Brandenberg and Stursberg 2021).

The reverse polar set was first introduced by Balas and Ivanescu (1964) and can be defined as follows.

Definition 3.

The reverse polar set of a set $S\subset {\mathbb{R}}^n$ is defined as $S^{-}≔\{d\in {\mathbb{R}}^n\hspace{0.17em}:\hspace{0.17em}{d}^{\top }x\le -1\hspace{0.17em}\forall x\in S\}$.

To simplify notation, we set ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)≔{(\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))-(x_{a(n)}^i,{\theta }_n^i))}^{-}$ for the reverse polar set of $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ shifted by $(x_{a(n)}^i,{\theta }_n^i)$. Using Remark 4, it can be reformulated.

Lemma 4.

The reverse polar set ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$ can be expressed as

$$\begin{array}{ll}\hfill & {\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)=\left\{({\gamma }_n,{\gamma }_{n0})\in {\mathbb{R}}^{d_{a(n)}}\times \mathbb{R}\hspace{0.17em}:\hspace{0.17em}\exists {\mu }_n\ge 0\hspace{0.17em}:\hspace{0.17em}\begin{array}{ll}\hfill & {\gamma }_{n0}\le 0,\hspace{0.17em}{\tilde{A}}_n^{\top }{\mu }_n-{\gamma }_{n0}{c}_n=0,\hspace{0.17em}-{\tilde{B}}_n^{\top }{\mu }_n-{\gamma }_n=0,\hfill \\ \hfill & \hspace{0.17em}{\tilde{d}}_n^{\top }{\mu }_n+{\gamma }_n^{\top }{x}_{a(n)}^i+{\gamma }_{n0}{\theta }_n^i\ge 1\hfill \end{array}\right\}.\hfill \end{array}$$

We provide a proof in Online Appendix EC.2.5.

Remark 8.

Even with the above reformulation of ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$, an explicit formulation is usually not readily available because of the existence quantor and because of ${\tilde{A}}_n,{\tilde{B}}_n$, and ${\tilde{d}}_n$ not being known.

Even though ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$ is not known explicitly, it proves useful in the derivation of Lagrangian cuts with favorable properties. This is because as we shall see, optimizing over ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$ is closely linked to solving normalized Lagrangian dual problems (11). So, in fact, our two aims for cut selection are intertwined and boil down to considering specific (bounded) normalizations of Problem (9).

In the next two subsections, we consider normalizations by norm constraints and by linear constraints. As we show, both approaches can be viewed from three different perspectives (a distance perspective, a projection perspective, and a reverse polar perspective) that help in understanding the properties of the related cuts. These perspectives have been analyzed in the literature before as stated in Table 2, but they have not been linked all together in a generalized framework and have not been applied specifically to Lagrangian cuts.

Table 2. Examination of the Normalized Cut-Generation Problems and Different Perspectives on It in the Literature

3.4. Normalization by Norm (Deep Lagrangian Cuts)

We consider the normalized Lagrangian dual (11) if some norm is used as the normalization function.

Definition 4.

Let $\Vert ·\Vert$ be some arbitrary norm. The Lagrangian cut (10) defined by the solution $({\pi }_n,{\pi }_{n0})$ to the normalized Lagrangian dual (11) with $g_n({\pi }_n,{\pi }_{n0})=\Vert {\pi }_n,{\pi }_{n0}\Vert$ is called a $\Vert ·\Vert$-deep Lagrangian cut. For ${\mathcal{\ell }}^p$-norms, we may also use the term ${\mathcal{\ell }}^p$-deep Lagrangian cuts.

If appropriate norms are used (e.g., ${\mathcal{\ell }}^1$ or ${\mathcal{\ell }}^{\infty }$), then the normalization can be expressed by linear constraints.

Deep cuts allow for three theoretical and geometrical interpretations (see Table 2). As the existing results from the literature can be applied to the multistage and Lagrangian setting in a straightforward way, for reasons of space, we do not provide proofs here.

1. Maximizing cut depth. Among all potential cuts, these cuts maximize the distance between the hyperplane associated with the cut and the incumbent $(x_{a(n)}^i,{\theta }_n^i)$ in the dual norm ${\Vert ·\Vert }_{*}$ of $\Vert ·\Vert$.

Lemma 5

(Based on Hosseini and Turner 2025). Let $\Vert ·\Vert$ be some norm and ${\Vert ·\Vert }_{*}$ be its dual norm. Further, let $d_n((x_{a(n)}^i,{\theta }_n^i);({\pi }_n,{\pi }_{n0}))$ denote the distance between the hyperplane defined by $({\pi }_n,{\pi }_{n0})$ and the point $(x_{a(n)}^i,{\theta }_n^i)\notin \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ measured in ${\Vert ·\Vert }_{*}$. Then, the optimal value ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)$ of Problem (11) with $g_n({\pi }_n,{\pi }_{n0})=\Vert {\pi }_n,{\pi }_{n0}\Vert$ equals

$$\underset{{\pi }_n,{\pi }_{n0}}{\max }\{d_n((x_{a(n)}^i,{\theta }_n^i);({\pi }_n,{\pi }_{n0}))\hspace{0.17em}:\hspace{0.17em}{\pi }_{n0}\ge 0\}.$$

• 2. Projection onto the epigraph. From a primal perspective, generating $\Vert ·\Vert$-deep cuts is in some sense equivalent to minimizing the distance in ${\Vert ·\Vert }_{*}$ between the incumbent $(x_{a(n)}^i,{\theta }_n^i)$ and the epigraph $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ (i.e., related to projecting $(x_{a(n)}^i,{\theta }_n^i)$ onto the epigraph).

Lemma 6

(Based on Cadoux 2010, Lemma 2.5). Let $\Vert ·\Vert$ be some norm and ${\Vert ·\Vert }_{*}$ be its dual norm. Then, the optimal value ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)$ of Problem (11) with $g_n({\pi }_n,{\pi }_{n0})=\Vert {\pi }_n,{\pi }_{n0}\Vert$ equals that of the projection problem

$$\begin{array}{l}\underset{x_{a(n)},{\theta }_n}{\min }\{{\Vert x_{a(n)}-x_{a(n)}^i,{\theta }_n-{\theta }_n^i\Vert }_{*}\hspace{0.17em}:\hspace{0.17em}(x_{a(n)},{\theta }_n)\in \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))\}.\hfill \end{array}$$(12)

• Lemma 6 implies ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)>0$ if $(x_{a(n)}^i,{\theta }_n^i)\notin \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ and ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)=0$ else, so it confirms Lemma 3. Therefore, as for the nonnormalized case, we have a unique flag for cases where no separating cut has to be constructed. However, in contrast to the nonnormalized case, the dual problem is bounded. We can also conclude from Lemma 6 that a deep cut supports $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$.

Corollary 1

(Based on Hosseini and Turner 2025, Proposition 3). Suppose $(x_{a(n)}^i,{\theta }_n^i)\notin \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$, and let $({\widehat{x}}_{a\left(n\right)},{\widehat{\theta }}_n)\in \text{epi}\left(\overline{\text{co}}\right({\underset{¯}{Q}}_n^{i+1}\left)\right)$ be a solution to (12). Then, any $\Vert ·\Vert$-deep cut separating $(x_{a(n)}^i,{\theta }_n^i)$ from $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ supports $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ at $({\widehat{x}}_{a\left(n\right)},{\widehat{\theta }}_n)$.

• 3. Minimizing a norm over the reverse polar set. Interestingly, deep cuts allow for another geometric interpretation that is related to the reverse polar set ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$.

Lemma 7

(Cadoux 2010, Lemma 2.9). Let $({\pi }_n^i,{\pi }_{n0}^i)$ be the coefficients of a $\Vert ·\Vert$-deep cut constructed at $(x_{a(n)}^i,{\theta }_n^i)\notin \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$. Then, there exists some $\alpha >0$ such that $-\alpha ({\pi }_n^i,{\pi }_{n0}^i)$ minimizes $\Vert ·\Vert$ over the reverse polar set ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$.

To illustrate these perspectives for different norms, we provide an example in Online Appendix EC.3. It highlights that whereas deep cuts can be unique, degenerate solutions are also possible if the optimization over ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$ in Lemma 7 does not have a unique solution. This degeneracy of the dual (11) may lead to selection of nondominant cuts (compare with challenge (b) in Section 1.1) and is only excluded for the ${\mathcal{\ell }}^2$-norm. Further, even if unique, deep cuts are not guaranteed to be facet defining. In fact, analyses in Cadoux (2010) and Hosseini and Turner (2025) show that they tend to be flat, at least when the optimality gap is still large. Finally, the projection problem (12) is also not guaranteed to have a unique solution for all but the ${\mathcal{\ell }}^2$-norm. However, if the solution is nonunique, the associated cut is unique and facet defining.

3.5. Linear Normalization (LN Lagrangian Cuts)

We consider the normalized Lagrangian dual (11) with a linear normalization function $g_n({\pi }_n,{\pi }_{n0})=u_n^{\top }{\pi }_n+u_{n0}{\pi }_{n0}$ defined by some coefficients $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \in \mathbb{R}$. Recall that we introduce the normalization constraint in (11) in order to transform Problem (9) to a bounded one. In contrast to the norm-based normalization from Section 3.4, a linear normalization does not necessarily imply boundedness. Hence, the choice of $(u_n,u_{n0})$ is crucial to ensure that an optimal solution exists. We further analyze this later in this section, but for now, we take the following assumption.

Assumption 3.

Given some $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \mathbb{R}$, the normalized Lagrangian dual (11) with $g_n({\pi }_n,{\pi }_{n0})=u_n^{\top }{\pi }_n+u_{n0}{\pi }_{n0}$ has a finite optimal value ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)>0$.

We can then define the associated type of Lagrangian cuts.

Definition 5.

Let $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \mathbb{R}$, and let the normalized Lagrangian dual (11) satisfy Assumption 3. Then, we refer to the associated cut (10) as a linear normalization (LN) Lagrangian cut.

Again, we can take three different perspectives on LN cuts.

1. Maximizing cut depth. As Hosseini and Turner (2025) show, LN cuts can be interpreted as maximizing the distance between the associated hyperplane and the incumbent in a more general type of distance function than a norm.

2. Projection on a line segment. Even though this geometric perspective has been discussed before for BD, most recently by Seo et al. (2022), we study it in more detail for the Lagrangian case as the formal description changes a bit. First, we exploit that for linear $g_n(·)$, the normalized Lagrangian dual (11) can be reformulated as an LP and then, be dualized with strong duality.

Lemma 8.

Let $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \mathbb{R}$. The normalized Lagrangian dual (11) with $g_n({\pi }_n,{\pi }_{n0})=u_n^{\top }{\pi }_n+u_{n0}{\pi }_{n0}$ can be formulated as an LP, and its dual is

$$\begin{array}{l}\underset{{\lambda }_n,z_n,{\eta }_n}{\min }\{{\eta }_n\hspace{0.17em}:\hspace{0.17em}({\lambda }_n,z_n)\in \text{conv}({\mathcal{W}}_n^{i+1}),\hspace{0.17em}{\eta }_n\ge 0,\hspace{0.17em}{u}_{n0}{\eta }_n\ge c_n^{\top }{\lambda }_n-{\theta }_n^i,\hspace{0.17em}{u}_n{\eta }_n=z_n-x_{a(n)}^i\}.\hfill \end{array}$$(13)

We provide a proof in Online Appendix EC.2.6. Problem (13) can be interpreted as finding the smallest scaling factor ${\eta }_n\ge 0$ such that starting from $(x_{a(n)}^i,{\theta }_n^i)$ along direction $(u_n,u_{n0})$, a point in $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ is reached. Again, for the optimal value, it follows ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)={\eta }_n^{*}>0$ if and only if $(x_{a(n)}^i,{\theta }_n^i)\notin \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$. Given the optimal value, the projection of $(x_{a(n)}^i,{\theta }_n^i)$ onto $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ along direction $(u_n,u_{n0})$ can be determined as

We then obtain the following result.

Corollary 2.

Let $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \mathbb{R}$, and let the normalized Lagrangian dual (11) satisfy Assumption 3. Furthermore, let $({\widehat{x}}_{a\left(n\right)},{\widehat{\theta }}_n)\in \text{epi}\left(\overline{\text{co}}\right({\underset{¯}{Q}}_n^{i+1}\left)\right)$ satisfy (14). Then, the associated LN cut supports $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ at $({\widehat{x}}_{a\left(n\right)},{\widehat{\theta }}_n)$.

• 3. Maximizing a linear function over the reverse polar set. Again, an interpretation with respect to ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$ is possible (see Cornuéjols and Lemaréchal 2006 and Brandenberg and Stursberg 2021 for details); LN cuts can be obtained by maximizing the linear objective function $u_n^{\top }{\gamma }_n+u_{n0}{\gamma }_{n0}$ over ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$:

  $$\underset{{\gamma }_n,{\gamma }_{n0}}{\max }\{u_n^{\top }{\gamma }_n+u_{n0}{\gamma }_{n0}\hspace{0.17em}:\hspace{0.17em}({\gamma }_n,{\gamma }_{n0})\in {\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)\}.$$(15)

  For the relation to the normalized Lagrangian dual (11), the following result holds. For a sketch of the proof, see Online Appendix EC.2.7.

Theorem 2

(Based on Brandenberg and Stursberg 2021). Let $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \mathbb{R}$, and let $(x_{a(n)}^i,{\theta }_n^i)\notin \text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$. Then, the following results hold.

1. If Problem (11) satisfies Assumption 3, then Problem (15) has a finite optimal value and vice versa.

2. The induced cuts for $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ are equivalent for both problems.

We illustrate the different perspectives on LN cuts again using an example in Online Appendix EC.3.

Geometrically, a favorable property of generating cuts based on ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$ is that (given that $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ is a full-dimensional polyhedron) each vertex $({\gamma }_n,{\gamma }_{n0})$ of ${\mathcal{R}}_n^{i+1}(x_{a(n)}^i,{\theta }_n^i)$ corresponds to the normal vector of a facet of $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ and vice versa (Brandenberg and Stursberg 2021). In order to identify such points, we can use the LP (15) given a choice of $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \mathbb{R}$ such that a finite optimum is attained. In fact, perspectives (2) and (3) make a sufficient condition for such a finite optimum readily available, and by that, they also allow us to conclude when Assumption 3 is satisfied. Whereas this result is known from convex analysis (Cornuéjols and Lemaréchal 2006, Brandenberg and Stursberg 2021), in Online Appendix EC.2.8, we provide an alternative proof based on perspective (2) and the specific Lagrangian setting.

Lemma 9.

Assumption 3 is satisfied if

$$(u_n,u_{n0})\in \text{cone}(\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))-(x_{a(n)}^i,{\theta }_n^i))\setminus \left\{0\right\},$$(16)

where $\text{cone}(S)$ denotes the conical hull of a set S.

Lemma 9 implies that also choosing $(u_n,u_{n0})$ from $\text{epi}({\underset{¯}{Q}}_n^{i+1})-(x_{a(n)}^i,{\theta }_n^i)$ or even from $\text{epi}(Q_n)-(x_{a(n)}^i,{\theta }_n^i)$ is sufficient. In other words, choosing reasonable coefficients $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \mathbb{R}$ boils down to finding a core point within one of these epigraphs. We discuss this in more detail in Section 3.6.

We now state some beneficial properties of LN cuts with respect to the aforementioned cut-quality criteria.

Theorem 3.

Let $(u_n,u_{n0})\in {\mathbb{R}}^{d_n}\times \mathbb{R}$. Consider the normalized Lagrangian dual problem (11) with $g_n({\pi }_n,{\pi }_{n0})=u_n^{\top }{\pi }_n+u_{n0}{\pi }_{n0}$.

1. For all $(u_n,u_{n0})\in \text{cone}(\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))-(x_{a(n)}^i,{\theta }_n^i))\setminus \{0\}$, the optimal point $({\pi }_n^{*},{\pi }_{n0}^{*})$ defines a supporting cut for $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$.

2. For all $(u_n,u_{n0})\in \text{cone}(\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))-(x_{a(n)}^i,{\theta }_n^i))\setminus \{0\}$, there exists an optimal extreme point $({\pi }_n^{*},{\pi }_{n0}^{*})$ such that the obtained cut is facet defining for $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$.

3. For all $(u_n,u_{n0})\in \text{relint}(\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))-(x_{a(n)}^i,{\theta }_n^i))$, any optimal point $({\pi }_n^{*},{\pi }_{n0}^{*})$ with ${\pi }_{n0}^{*}>0$ defines a Pareto-optimal cut for $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ on $\text{conv}(Z_{a(n)})$.

Case (i) of Theorem 3 directly follows from Lemma 9 and Corollary 2. Case (ii) of Theorem 3 follows from Brandenberg and Stursberg (2021, theorem 3.3), and case (iii) of Theorem 3 follows from Stursberg (2019, theorem 3.43). Note that the results from the literature require us to choose $(u_n,u_{n0})$ from the relative interior of the epigraph restricted to $\text{conv}(X_{a(n)})$. However, under Assumption 1, if we choose $Z_n=X_{a(n)}$ or $Z_n=\text{conv}(X_{a(n)})$, in our case, the considered epigraphs are always restricted to this set.

Considering case (ii) of Theorem 3, the only case in which no facet-defining cut is obtained occurs if the optimal solution to the normalized Lagrangian dual (11) is not unique. However, this is only the case for a small subset of choices $(u_n,u_{n0})$. With respect to case (iii) of Theorem 3, we should emphasize again that Pareto optimality for $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ does not necessarily imply Pareto optimality for $\text{epi}({\underset{¯}{Q}}_n^{i+1})$.

3.6. Identifying Core Points

A key challenge of generating LN cuts with favorable properties is to choose $(u_n,u_{n0})$ appropriately (i.e., according to Lemma 9 or Theorem 3). In the literature, it is often proposed to evaluate feasible points in the true objective function to obtain core points $({\check{x}}_{a(n)}^i,{\check{\theta }}_n^i)$ and by that, coefficients $(u_n,u_{n0})$ that best guide the solution process. Whereas this approach is straightforward for BD or the two-stage stochastic case (Stursberg 2019), in the multistage case, evaluating $Q_n(·)$ exactly is computationally prohibitive in general. Therefore, we propose different heuristics based on using function ${\underset{¯}{Q}}_n^{i+1}(·)$.

• Mid. We assume that $x_{a(n)}$ has to satisfy box constraints and then, set ${\check{x}}_{a(n)}^i$ to the component-wise midpoint. For instance, if $X_{a(n)}={\{0,1\}}^{d_{a(n)}}$, then we set ${\check{x}}_{a(n)}^i=\text{mid}(\text{conv}(X_{a(n)}))={(\frac{1}{2},\dots ,\frac{1}{2})}^{\top }$. We then set ${\check{\theta }}_n^i={\underset{¯}{Q}}_n^{i+1}({\check{x}}_{a(n)}^i)$, even if this does not necessarily satisfy the relative interior requirement in Theorem 3. The idea of this approach is that incentivizing the LN cuts to support $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ in the interior of the (convexified) state space may be useful to avoid the degeneracy issues discussed in challenge (b) in Section 1.1.

• Eps. For some $\epsilon >0$, we set ${\check{x}}_{a(n)}^i$ to an $\epsilon$-perturbation of $x_{a(n)}^i$ into the interior of $\text{conv}(X_{a(n)})$. For instance, if $X_{a(n)}={\{0,1\}}^{d_{a(n)}}$, then we set ${\check{x}}_{a(n),k}^i=1-\epsilon$ for each k such that $x_{a(n),k}^i=1$ and ${\check{x}}_{a(n),k}^i=\epsilon$ for each k such that $x_{a(n),k}^i=0$. Again, we set ${\check{\theta }}_n^i={\underset{¯}{Q}}_n^{i+1}({\check{x}}_{a(n)}^i)$. This approach is inspired by the perturbation strategy described in Sherali and Lunday (2013) and Seo et al. (2022). It is particularly suited to SDDiP because it avoids the possible degeneracy issues related to generating cuts at extreme points of the state space; see Section 1.1. It should thus favor the generation of facet-defining cuts, which for sufficiently small $\epsilon$, should still be supporting $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ at $(x_{a(n)}^i,\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1})(x_{a(n)}^i))$.

• Relint. Based on Lemma EC.1 in Online Appendix EC.2.1, we may solve feasibility problem (7) but replace parameter $x_{a(n)}^i$ with variable $x_{a(n)}$, and then, we set $({\check{x}}_{a(n)}^i,{\check{\theta }}_n^i)$ to any feasible point obtained. However, in the light of Theorem 3, to incentivize finding a point in $\text{relint}(\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1})))$, we additionally introduce tightening variables and then, maximize the number of inactive constraints. For technicals details, we refer to our implementation and to Conforti and Wolsey (2019, section 5.1), where a similar strategy is used.

• Conv. We note that any convex combination of two (or more) feasible points in Subproblem (6) is always contained in $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$. One such point is readily available with ${\rho }_1≔(x_{a(n)}^i,{\underset{¯}{Q}}_n^{i+1}(x_{a(n)}^i))$. For the case $X_{a(n)}={\{0,1\}}^{d_{a(n)}}$, an intuitive strategy to obtain a second point ${\rho }_2$ is to consider the diagonal counterpart of $x_{a(n)}^i$ (i.e., swapping zero and one for all components) and the associated function value for ${\underset{¯}{Q}}_n^{i+1}(·)$. As this point, ${\rho }_2$ is feasible under Assumption 1, and we can compute a core point as $({\check{x}}_{a(n)}^i,{\check{\theta }}_n^i)=\alpha {\rho }_1+(1-\alpha ){\rho }_2$ for some predefined $\alpha \in (0,1)$.

The Mid, Eps, and Conv heuristics are visualized in Figure 1. Note that for the given example and Eps, the generated cut is facet defining and corresponds to the best-possible tight cut that can be obtained using the standard cut generation from SDDiP. However, in contrast to the latter approach, the cut-generation problem is not degenerate, so there is no risk of generating tight but much weaker cuts. Of course, in the general multistage case, multiple different facets may exist, and Eps is not guaranteed to choose the best one.

Despite their straightforwardness, these heuristics come with some notable challenges. The first two heuristics yield candidates satisfying ${\check{x}}_{a(n)}^i\in \text{conv}(X_{a(n)})$. However, this choice is not necessarily feasible if integrality constraints are present (especially if $Z_{a(n)}=X_{a(n)}$), leading to ${\underset{¯}{Q}}_n^{i+1}({\check{x}}_{a(n)}^i)=+\infty$ and preventing a reasonable choice for $(u_n,u_{n0})$. Instead of evaluating ${\underset{¯}{Q}}_n^{i+1}(·)$, in such a case, we may revert to the value of the associated LP relaxation (LP-value), the current state component ${\theta }_n^i$ (epi-state), or the objective value ${\underset{¯}{Q}}_n^{i+1}(x_{a(n)}^i)$ at the incumbent (primal-obj). However, none of these approaches are guaranteed to yield core points. The heuristic Conv does not face this exact issue as it always considers convex combinations of objective values of feasible points. On the other hand, to be able to identify a good core point using this heuristic, it is first required to find a reasonable second feasible solution ${\rho }_2$. This can be challenging in practice if relatively complete recourse is not satisfied or if $X_{a(n)}\ne {\{0,1\}}^{d_{a(n)}}$. We provide a visualization of these challenges and improvement strategies in Online Appendix EC.4.

For these reasons, even with the proposed heuristics, we may end up with a normalized Lagrangian dual problem (11) that is unbounded. This is clearly unintended and should be detected in practice. We discuss strategies to do this in Online Appendix EC.5. Whenever (potential) unboundedness is detected, we may take some special countermeasures, such as artificially bounding the normalized Lagrangian dual problem (11) or resorting to strengthened Benders cuts instead of generating LN cuts.

4. Convergence Aspects

Given purely binary-state variables (i.e., $x_n\in {\{0,1\}}^{d_n}$ for all $n\in \mathcal{N}$), the classical Lagrangian cuts from Zou et al. (2019) are valid, tight for ${\underset{¯}{Q}}_n^{i+1}(·)$ at the incumbent $x_{a(n)}^i$, and finite (see Online Appendix EC.1). Therefore, they are sufficient to ensure (almost sure) finite convergence of SDDiP. Whereas they are always supporting $\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1})(·)$ according to Corollaries 1 and 2, deep and LN Lagrangian cuts are not necessarily tight at the incumbent, so these convergence results are not directly applicable. However, in this section, we show that under certain assumptions, the finite convergence results of SDDiP can still be established. For reasons of space, we do not present a complete, self-contained convergence proof but try to convey the key ideas.

4.1. Convergence for LN Cuts

For LN cuts, recall that according to Theorem 3, they are almost surely facet defining for $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ (given that $(u_n,u_{n0})$ satisfies the conditions of case (ii) of Theorem 3). Additionally, $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ is a convex polyhedron with finitely many facets (see Remark 4 and the formulation in Online Appendix EC.5). Therefore, if in all iterations (or at least at finite intervals of iterations), facet-defining LN cuts are generated, at most finitely many different cuts can be generated before $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ is recovered completely. For the binary state space, this implies a tight approximation of ${\underset{¯}{Q}}_n^{i+1}(·)$ at all $x_{a(n)}\in {\{0,1\}}^{d_{a(n)}}$. This is sufficient to obtain the known convergence results for SDDiP. For more technical details, we refer to Füllner (2024, chapter D.4).

4.2. Convergence for Deep Cuts

For deep Lagrangian cuts, proving convergence is a bit more intricate. The main idea, which is visualized in Online Appendix EC.2.9, is the following. Even if these cuts are not guaranteed to be tight at $x_{a(n)}^i$ or facet defining, tight cuts will eventually be generated after finitely many steps. To show this, we use two key concepts. First, we consider subproblems with fixed approximations ${\mathrm{\Psi }}_m$ for all $m\in \mathcal{C}(n)$ and all $n\in \overline{\mathcal{N}}$. In this case, the set $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))$ remains constant over the iterations, so we drop its iteration index for simplicity. This assumption is satisfied for all leaf nodes in $\mathcal{N}$. All of the results that we derive for this case can then be extended to the whole scenario tree $\mathcal{T}$ using induction. Second, we consider a subsequence of incumbents with fixed first component ${\overline{x}}_{a(n)}$. If $X_n={\{0,1\}}^{d_n}$ for all $n\in \mathcal{N}$, then only finitely many different values can be attained by the state variables. Therefore, all of the results that we derive in this section can be extended to the whole state space by finite repetition.

For notational convenience, we set ${\overline{\theta }}_n≔\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1})({\overline{x}}_{a(n)})$. We first obtain the following auxiliary result.

Lemma 10.

For any $n\in \overline{\mathcal{N}}$, consider a subsequence of incumbents ${({\overline{x}}_{a(n)},{\theta }_n^{i_{\nu }})}_{\nu \in \mathbb{N}}$ with fixed first component ${\overline{x}}_{a(n)}\in X_{a(n)}$ for all $\nu \in \mathbb{N}$. Then, the point $({\overline{x}}_{a(n)},{\overline{\theta }}_n)$ and any sequence ${({\widehat{x}}_{a(n)}^{i_{\nu }},{\widehat{\theta }}_n^{i_{\nu }})}_{\nu \in \mathbb{N}}$ of solutions of the projection problems (12) satisfy ${\lim }_{\nu \to \infty }({\widehat{x}}_{a(n)}^{i_{\nu }},{\widehat{\theta }}_n^{i_{\nu }})=({\overline{x}}_{a(n)},{\overline{\theta }}_n)$.

We provide a proof in Online Appendix EC.2.10. Now, let K with $|K|\in \mathbb{N}$ denote the finite set of facets F of $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n))$. Furthermore, let $\overline{K}\subseteq K$ denote the subset of facets such that $({\overline{x}}_{a(n)},{\overline{\theta }}_n)\in F_k$ for all $k\in \overline{K}$ and $({\overline{x}}_{a(n)},{\overline{\theta }}_n)\notin F_k$ for any $k\notin \overline{K}$, and analogously, let ${\widehat{K}}^{\nu }\subseteq K$ denote the subset of facets in which $({\widehat{x}}_{a(n)}^{i_{\nu }},{\widehat{\theta }}_n^{i_{\nu }})$ is contained. By definition of the points, $\overline{K}$ and ${\widehat{K}}^{\nu }$ are nonempty. Intuitively, from Lemma 10, we can conclude that for sufficiently large $\nu$, the points $({\widehat{x}}_{a(n)}^{i_{\nu }},{\widehat{\theta }}_n^{i_{\nu }})$ and $({\overline{x}}_{a(n)},{\overline{\theta }}_n)$ are located on at least one joint facet and that $({\widehat{x}}_{a(n)}^{i_{\nu }},{\widehat{\theta }}_n^{i_{\nu }})$ can only be a vertex of $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n))$ if both points coincide. More precisely, we obtain the following result. For a detailed proof, we refer again to Füllner (2024, chapter D.4).

Lemma 11.

There exists some $\widehat{\nu }\in \mathbb{N}$ such that for all $\nu \ge \widehat{\nu }$, we have ${\widehat{K}}^{\nu }\subseteq \overline{K}$. Moreover, there exists some $\overline{\nu }\ge \widehat{\nu },\overline{\nu }\in \mathbb{N},$ such that for all $\nu \ge \overline{\nu }$, the point $({\widehat{x}}_{a(n)}^{i_{\nu }},{\widehat{\theta }}_n^{i_{\nu }})$ either is equal to $({\overline{x}}_{a(n)},{\overline{\theta }}_n)$ or is not a vertex of $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n))$.

Moreover, as also shown in Füllner (2024, chapter D.4), cuts supporting $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n))$ at a point in the relative interior of a face $\mathfrak{F}$ (or facet F) are also tight approximators at all other points at the same face (or facet). Therefore, because of Lemma 11, cuts supporting $({\widehat{x}}_{a(n)}^{i_{\nu }},{\widehat{\theta }}_n^{i_{\nu }})$ eventually will also support $({\overline{x}}_{a(n)},{\overline{\theta }}_n)$. With these auxiliary results, we are able to state the main tightness result, which we prove in Online Appendix EC.2.11.

Lemma 12.

For any $\overline{\mathcal{N}}$, consider Subproblem (6) with fixed approximations ${\mathrm{\Psi }}_m$ for all $m\in \mathcal{C}(n)$. Moreover, consider a subsequence of trial points ${({\overline{x}}_{a(n)},{\theta }_n^{i_{\nu }})}_{\nu \in \mathbb{N}}$ with fixed first component ${\overline{x}}_{a(n)}\in X_{a(n)}$ for all $\nu \in \mathbb{N}$. Then, there exists some $\check{\nu }\in \mathbb{N}$ such that for all $\nu \ge \check{\nu }$, a deep Lagrangian cut as defined in Definition 4 supports $\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n))$ at $({\overline{x}}_{a(n)},{\overline{\theta }}_n)$.

We should emphasize that this worst-case reasoning of deep cuts eventually coinciding with classical Lagrangian cuts is solely used to show that (almost sure) finite convergence of SDDiP is still assured. In practice, the idea is that deep (as well as LN) Lagrangian cuts lead to stronger approximations outside of $x_{a(n)}^i$ and thus, improve the computational performance.

5. Computational Experiments

We report results for a computational study of SDDiP incorporating the proposed cut-generation framework (Algorithm 1). For a detailed description of the SDDiP algorithm, we refer to Zou et al. (2019). We run SDDiP for two classes of MS-MILPs from the literature.

• Capacitated lot-sizing problem (CLSP). A capacitated lot-sizing problem described in Trigeiro et al. (1989) with stagewise independent uncertain demand for each product. This problem is also considered in Ahmed et al. (2022) and identified to be challenging for exact decomposition methods, like SDDiP. We consider instances with 3 or 10 state variables; 4, 10, or 16 stages, and 20 realizations of the uncertain demand at each stage.

• Capacitated facility location problem (CFLP). A capacitated facility location problem with stagewise independent uncertain facility disruptions (Seranilla and Löhndorf 2024). We consider an instance with 16 state variables, 100 stages, and 20 realization of the uncertainty at each stage.

With respect to the state space, it is important to mention that CLSP contains continuous state variables. For instances with three state variables, in order to ensure exactness of SDDiP, we thus apply a binary approximation of the state space using a discretization precision of 1.0; see Zou et al. (2019). This means that the modified MS-MILP, which is tackled by SDDiP, contains 30 state variables. We refer to this problem as CLSP-Bin. For instances with 10 state variables, using a binary approximation would produce state dimensions that are computationally intractable for SDDiP as its complexity grows exponentially in the dimension of the state space. Therefore, in these cases, we apply SDDiP to CLSP without a binary approximation and thus, without theoretical guarantees to close the optimality gap. In contrast, CFLP only contains binary state variables from the outset, so no binary approximation is required.

With respect to the uncertainty, for CLSP-Bin instances, we use the exact same scenarios as in Ahmed et al. (2022). For the larger instances (CLSP) that are not covered in that paper, we generate new scenarios using the same methodology. For CFLP, we use scenarios provided to us by Seranilla and Löhndorf (2024).

To assess the quality of the cuts obtained using our framework, we also run comparative tests using established cut-generation techniques in SDDiP. More precisely, we consider the following types of cuts.

• B, SB, and L. Classical Benders cuts, strengthened Benders cuts (Zou et al. 2019), and Lagrangian cuts (Zou et al. 2019) (see Online Appendix EC.1), respectively, using a single-cut or multicut approach.

• ${\mathcal{\ell }}^1$ and ${\mathcal{\ell }}^{\infty }$. Deep Lagrangian cuts from Definition 4 for the ${\mathcal{\ell }}^1$-norm and the ${\mathcal{\ell }}^{\infty }$-norm, respectively.

• LN. LN Lagrangian cuts from Definition 5 using the Mid, Eps, Relint, and Conv heuristics from Section 3.6 to identify core points and to determine the normalization coefficients.

The test cases and cut strategies are summarized in Table 3. By construction, all deep and LN cuts require using a multicut approach. For Eps, we use a perturbation of the incumbent toward $\text{mid}(\text{conv}(X_{a(n)}))$ by ${10}^{-2}$ in each component. Conv is only considered for problem CFLP with binary state space as then, ${\rho }_2$ can be computed in a straightforward way; see Section 3.6. We use values of 0.5, 0.75, 0.9, and 0.99 for $\alpha$. Additionally, for problem CFLP, the risk of ${\underset{¯}{Q}}_n^{i+1}({\check{x}}_{a(n)}^i)=+\infty$ is prevalent in the search for core points because of the binary nature of the state space. Therefore, we use the objective value of the LP relaxation instead (LP-value), and we apply the strategy primal-obj to lift the $\theta$-component of the core-point candidate as described in Section 3.6.

Table 3. Overview of Test Problems and Configuration of Experiments

5.1. Implementation Details

SDDiP and all cut-generation approaches are implemented in Julia-1.9.3 (Bezanson et al. 2017) based on the existing packages SDDP.jl (Dowson and Kapelevich 2021) and JuMP.jl (Dunning et al. 2017). The code is available in GitHub (Füllner et al. 2026).

SDDiP is terminated after a predefined time limit or if the lower bounds do not improve considerably for 20 iterations. In each forward pass, one scenario path is randomly sampled. After termination of SDDiP, an in-sample Monte Carlo simulation with 1,000 replications is conducted on the finite scenario tree to compute a statistical upper bound (UB) for the obtained policy.

For LN cuts, the choice of normalization coefficients $(u_n,u_{n0})$ is crucial for the cut quality but also, to achieve a bounded dual problem. For all of the heuristics presented in Section 3.6, we may obtain an unbounded problem (11) if $(u_n,u_{n0})\notin \text{cone}(\text{epi}(\overline{\text{co}}({\underset{¯}{Q}}_n^{i+1}))-(x_{a(n)}^i,{\theta }_n^i))\setminus \{0\}$; see Lemma 9. To detect a potentially insufficient choice of $(u_n,u_{n0})$ in advance, we use a (conservative) restriction of the problem (Figure EC.9 in Online Appendix EC.7) and use its infeasibility as a flag for possible unboundedness. This strategy is explained in more detail in Online Appendix EC.5. Given such a flag, we only generate an SB cut instead of an LN cut.

The Lagrangian dual problems in SDDiP are solved with a level bundle method (Lemaréchal et al. 1995), which can be considered as an enhancement of the Kelley method (Kelley 1960) (Online Appendix EC.6). Additionally, we bound the dual objective by ${\underset{¯}{Q}}_n^{i+1}(x_{a(n)}^i)$. We run the method for a maximum of 1,000 (for CLSP) or 100 (for CFLP) iterations and an optimality tolerance of ${10}^{-4}$ (for CLSP) or ${10}^{-2}$ (for CFLP). The multipliers ${\pi }_n$ are initialized with a vector of zeros, and ${\pi }_{n0}$ is initialized with one. If the level bundle method reports infeasibilities in the quadratic auxiliary problem, we proceed with a standard Kelley step instead. In the event of other numerical issues, a valid cut is constructed with the current values of the multipliers.

Additionally, in order to avoid numerical issues from the generated cuts, we apply some countermeasures.

• We only add a new cut if the previous incumbent $(x_{a(n)}^i,{\theta }_n^i)$ is cut away by at least ${10}^{-4}$.

• We do not add a cut if $(u_n,u_{n0})\approx \mathbf{0}$, if ${\pi }_{n0}\approx 0$, or if ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)\approx 0$.

For CLSP-Bin and approach ${\mathcal{\ell }}^{\infty }$, some preliminary tests indicated a proneness to degeneracy of the dual problem and thus, cuts of bad quality. Therefore, for this case, we perform an additional optimization step where we minimize the ${\mathcal{\ell }}^1$-norm of ${\pi }_n$ among all optimal dual solutions (to which we restrict by fixing the dual objective value to ${\widehat{v}}_n^{ND,i+1}(x_{a(n)}^i,{\theta }_n^i)$). This is indicated by MNC (minimal norm choice) in Table 3. This approach was previously proposed in SDDP.jl (Dowson and Kapelevich 2021). It is not applied for CLSP and CFLP.

All occurring LP, MILP, and quadratic subproblems are solved using Gurobi 9.0.3 with an optimality tolerance of ${10}^{-4}$ and a time limit of 300 seconds. All tests are run on a Windows machine with 64 GB RAM and an Intel Core i7-7700K processor (4.2 GHz).

5.2. Results

For our first batch of experiments, we apply SDDiP with only one type of cut per run. The results are presented in Tables 4–6 and Online Appendix EC.7. The columns in Tables 4–6 contain the number of stages, the used cut-generation approach (Method), the best lower bound obtained (Best LB), a statistical upper bound computed after termination of SDDiP (Stat. UB; we report the upper limit of the computed confidence interval), the optimality gap (Gap), the number of iterations (# Iter.), the time in seconds (Time), the average time per iteration (Time/It), the average number of iterations (Avg. Lag-It) and the time (Time/Lag) required in the level bundle method to solve the Lagrangian dual, the average number of generated cuts per node (Cuts/Stage), and for LN cuts, the share of cases with detected potential unboundedness in the core-point identification process (Unb). In some cases, the simulation did not yield an upper bound estimate because of numerical issues.

Table 4. SDDiP Results for CLSP-Bin

Table 5. SDDiP Results for CLSP with 10 State Variables

Table 6. SDDiP Results for CFLP

5.2.1. CLSP-Bin.

In addition to Table 4, the development of the lower bounds obtained within SDDiP is illustrated in Figure 2 for CLSP-Bin. For $T=4$, we can see that using deep or LN Lagrangian cuts, the LBs converge to the true optimal value. They clearly outperform the LBs obtained using approaches B, SB, and L, for which the bounds stall very fast and far away from optimality. In particular, despite the theoretical guarantees of the cuts, for L, the LBs improve very slowly for reasons outlined in Section 1.

For larger T, we observe that LN cuts perform significantly better than deep cuts, which tend to be very flat; see Online Appendix EC.7. We also detect no unboundedness of the dual (11), so the core-point heuristics seem to work reasonably well. Again, using LN cuts, the LBs from B, SB, and L are outperformed. However, compared with SB, the difference gets smaller the more stages T are considered. This is mostly for computational reasons. As stated in Table 4, many fewer iterations are finished within the time limit because of the high computational cost of solving Lagrangian dual problems for each stage and iteration (see columns Avg. Lag-It and Time/Lag in Table 4 and also, see Online Appendix EC.7). Compared with single-cut approaches, the number of cuts per subproblem also grows more quickly, further slowing down the solution process.

The LBs give an indication of the direct effect that different types of cuts have on the outer approximation quality and the convergence behavior of SDDiP. However, in practice, what we are most interested in is the quality of the feasible policies obtained from running SDDiP (i.e., the statistical upper bounds computed after termination of SDDiP). In this regard, whereas using deep or LN cuts may considerably improve the solution quality compared with L, this is often not true compared with SB, where good feasible solutions seem to be identified quickly, and many more iterations are finished within the time limit.

5.2.2. CLSP.

For experiments on CLSP, we do not apply a binary approximation, and thus, we cannot expect the optimality gap to be closed. However, in return, iterations should take considerably less time.

In fact, we observe that the average iteration time is reduced considerably (by 60%–75%) compared with CLSP-Bin. For this reason and because the considered state space is lower dimensional, even without theoretical convergence guarantees, better optimality gaps are obtained for CLSP in the same time.

Again, the results show that deep and LN Lagrangian cuts manage to achieve better LBs and optimality gaps than conventional cuts as presented in Figure 3 and Table 5. Overall, deep cuts perform better as more iterations are finished in the same time. The results show that our proposed techniques for generating Lagrangian cuts may be helpful to improve the convergence behavior of SDDiP, even if no binary approximation is applied. However, as for the previous experiments, we cannot conclude that the improvement in lower bounds necessarily leads to an improvement of the in-sample performance of the obtained policy.

5.2.3. CFLP.

For our experiments on CFLP, with binary state variables from the outset, the results are presented in Figure 4 and Table 6. The results differ considerably from those for the previous experiments.

First, we observe that all types of Lagrangian cuts lead to better UBs and thus, more feasible policies than SB. Interestingly, this is the case, although SB cuts are tight for this test problem. In reverse, deep cuts do not yield better LBs than B, SB, or L. This is mainly because deep cuts are completely flat, and thus, they lead to an unfavorable exploration of the state space (this also explains the small number of required iterations to solve the dual (11) to optimality). Concrete numbers to back these observations are provided in Online Appendix EC.7.

Second, for LN cuts, Eps yields tight cuts, but they seem to relate to very steep facets of the epigraph (again, see Online Appendix EC.7), so the obtained LBs are pretty bad. In contrast, Mid and Relint do not yield tight cuts on average. Additionally, for both, we observe the infeasibility challenges discussed in Section 3.6. They lead to (potentially) insufficient choices of $(u_n,u_{n0})$, and thus, potential unboundedness of Problem (11) for more than 80% of the cases so that SB cuts are generated instead. Still, generating alternative cuts in only 15%–20% of the cases improves the LBs and UBs by a lot.

In general, LN cuts that use core points with ${\check{x}}_{a(n)}^i$ relatively close to $\text{mid}(\text{conv}(X_{a(n)}))$ in early iterations (i.e., Mid, Conv(50), and Conv(75)) yield extremely good lower bounds in very few iterations and almost manage to close the optimality gap. According to Online Appendix EC.7, they seem to achieve a good trade-off between cut tightness at the incumbent and approximation quality everywhere else. Overall, the Conv heuristic leads to the best LBs and UBs without the infeasibility challenges that come along with Mid or Relint. The only drawback is that the average time spent for one iteration is extremely high.

5.3. Combination with Strengthened Benders Cuts

Using SDDiP with only Lagrangian cuts becomes extremely slow for large problems. Therefore, as already proposed in Zou et al. (2019), it is reasonable to combine them with cheaper cuts. In addition, the experiments for CLSP-Bin and CLSP indicate that using SB, it is often possible to quickly identify good feasible solutions but that the LBs may be too loose to get a certificate for optimality, whereas for our proposed Lagrangian cuts, it is the opposite. For these reasons, we consider a second batch of experiments where we start with only SB for the first 20 iterations to get a quick bound improvement and then, generate both SB cuts and Lagrangian cuts in each iteration. The results are presented in Figure 4 and Tables 4–6.

The LBs heavily improve for L but are not affected too much for deep or LN cuts. They are still better, however, than the ones obtained using only SB. As the simulated UBs are better than in the previous setting for all types of cuts, we conclude that a combination of SB cuts and Lagrangian cuts may be helpful to identify good policies early while not compromising good LBs and convergence guarantees.

5.4. Discussion and Potential Improvements

Overall, our results show significant improvements of the obtained LBs and UBs using the proposed cut-generation framework compared with standard Lagrangian cuts from Zou et al. (2019). In the light of challenges (b) and (c) presented in Section 1, this confirms that either better tight cuts or nontight cuts that improve the exploration of the state space are generated. For problem CFLP, the proposed LN cuts manage to approximately close the optimality gap, which is by far not the case for conventional cuts.

As advice for practitioners, it seems that deep cuts, which are often relatively flat, are most valuable for problems with large (initially continuous) state spaces as they are less costly to generate than LN cuts, which usually yield better LBs. LN cuts are more prone to numerical issues and require more tuning (e.g., with regard to the core-point heuristics). Among deep cuts, using the ${\mathcal{\ell }}^1$-norm seems most promising.

For (initially) continuous state spaces, using the heuristics Mid, Relint, and Eps yields very good LBs, with Mid looking like the most reliable and easy to implement approach. Also, almost no cases of potential unboundedness of the dual (11) are detected. This confirms that these heuristics work reasonably well for these types of problems. For problems with binary state spaces, however, infeasibility issues significantly aggravate the identification of (good) core points. In contrast, the Conv heuristic is tailor made for these problems. Under relatively complete recourse, it only requires one additional evaluation of ${\underset{¯}{Q}}_n^{i+1}(·)$ to compute a proven core point. As shown for CFLP, this strategy can be costly but can lead to very good results.

Whereas deep and LN variants outperform classical Lagrangian cuts, this is not necessarily true for SB cuts. For CLSP-Bin and CLSP, SB cuts yield worse LBs but better feasible policies at termination. Given their cheap computation, this means that SB cuts should remain the first choice for many practical problems, especially if SDDiP has to be stopped prematurely. If they do not close the optimality gap sufficiently, they may be combined with our proposed cuts to find a good balance between improving LBs and UBs.

We should acknowledge that our proposed framework does not necessarily lead to a reduced time spent on generating Lagrangian cuts. It may even be increased because of the requirement of using a multicut approach. Hence, challenge (a) from Section 1 remains. In the future, this could be addressed in several ways. First, cut-selection schemes, where the number of cuts is kept at a minimum, could be explored. Second, parallelization, warm-starting, or acceleration techniques for the Lagrangian dual could be explored. Third, a smart restriction of the dual feasible space (see Chen and Luedtke 2022) may help to reduce the computational effort for solving Lagrangian dual problems while not compromising cut quality by too much. Finally, it might be of interest to investigate potential extensions of the proposed framework to a single-cut approach.

6. Conclusion

In this article, we show how Lagrangian cuts can be generated by solving special normalized Lagrangian dual problems when solving MS-MILPs. We prove that different normalizations lead to cuts with different properties. This allows for a lot of flexibility and notably, extends the toolbox for cut generation in SDDiP.

We provide computational results for experiments on two known test problems. They show that lower bounds and also, sometimes upper bounds obtained in SDDiP can be vastly improved. However, other known drawbacks of SDDiP persist, such as excessive computational effort, slow convergence, and inability to close the optimality gap for large problem instances. Therefore, more research on cut selection, cut combination schemes, or dual space restrictions is required to allow for an extensive application in practice.

Finally, it should be possible to extend the proposed framework to cases where Lagrangian duals are solved to generate nonconvex cuts; see, for instance, Füllner and Rebennack (2022) and Füllner et al. (2024).