\begin{equation}gin{document}

\title{{\TheTitle}

\begin{equation}gin{abstract} In this paper we investigate  the dual of a Multistage
Stochastic Linear Program (MSLP) to study two related questions for this class of problems. The first of these questions is the study of the
optimal value of the problem as a function of the involved  parameters.

For this sensitivity analysis problem, we provide formulas for the derivatives of the value function with respect to the parameters and illustrate their application
on an inventory problem. Since these formulas involve optimal dual solutions,
we need an algorithm that computes such solutions to use them, i.e., we need to solve
the dual problem.

In this context, the second question we address is the study of solution
methods for the dual problem. Writing Dynamic Programming equations for the dual, we
can use an SDDP type method, called Dual SDDP, which solves these Dynamic Programming equations
computing a sequence of nonincreasing deterministic upper bounds on the optimal value
of the problem. However, applying
this method will only be possible if the Relatively Complete Recourse (RCR) holds for the dual.
Since the RCR assumption may fail to hold (even for simple problems),
we design two variants of Dual SDDP, namely Dual SDDP with penalizations and Dual SDDP with
feasibility cuts, that converge to the optimal value of the dual (and therefore primal when
there is no duality gap) problem under mild assumptions. We also show that optimal dual
solutions can be obtained computing dual solutions of the subproblems solved when applying
Primal SDDP to the original primal MSLP.

The study of this second question allows us to take a fresh look at the
class of MSLP with interstage dependent cost coefficients. Indeed, for this
class of problems, cost-to-go functions are non-convex and solution methods were so far
using SDDP for a Markov chain approximation of the cost coefficients process.
For these problems, we propose to apply Dual SDDP with penalizations to the
cost-to-go functions of the dual which are concave.
This algorithm converges to the optimal value of the problem.

Finally, as a proof of concept of the tools developed, we present the results
of numerical experiments computing the sensitivity of the optimal value of an inventory
problem as a function of parameters of the demand process
and compare Primal and Dual SDDP on the inventory and a hydro-thermal planning problems.
\varepsilonnd{abstract}

\begin{equation}gin{keywords}
Stochastic optimization, Sensitivity analysis, SDDP, Dual SDDP, Relatively complete recourse.
\varepsilonnd{keywords}

\begin{equation}gin{AMS}
 90C15, 90C90, 90C30
\varepsilonnd{AMS}

\mathcal{S}ection{Introduction}\label{sec-intr}

Duality plays a key role in optimization. For generic optimization problems,
weak duality allows  to bound the optimal value.
Dual information is also used in many optimization algorithms such as
Uzawa algorithm \cite{uzawabook}, primal-dual projected gradient \cite{nedic} or Stochastic Dual Dynamic Programming (SDDP) \cite{pereira}. Moreover, for several classes of optimization problems, the dual is easier to solve than the primal problem, for instance
when it is amenable to decomposition techniques such as price decomposition \cite{bertsekas}.
Even when there is a duality gap between the primal and dual optimal values, solving the dual
already gives a bound on the optimal value, as mentioned earlier.
Duality is also a fundamental tool in the reformulation of Robust Optimization problems, see for instance \cite{nem2}.
Finally, derivatives of the value function of classes of optimization problems can be
related to optimal dual solutions, see \cite{Bonnans2000}, \cite{Rockafellar}
and more recently \cite{guiguessiopt2016, guiguesinexactsddp, guiguesinexactsmd} for
the characterization of subdifferentials, subgradients, and
${\rm Var}epsilon$-subgradients of value functions of convex optimization problems.

For stochastic control problems, stochastic Lagrange multipliers
were already used in \cite{kushner69, kushner65, kushner72}.
In the context of multistage stochastic programs, duality was studied in
\cite{rockduality99, higlesen2006}, see also \cite{shadenrbook} for a review.
More recently, the sensitivity analysis of multistage stochastic programs was discussed in \cite{bonnan}  and \cite{woza}.
In \cite{bonnan} the authors study the sensitivity with respect to parameters driving the considered  price model. The corresponding parameters are in the objective function and the analysis of  the estimate of marginal price is  based on  Danskin's theorem with the SDDP
method used for the numerical calculations.
In  \cite{woza}, the authors use the Envelope Theorem for the sensitivity analysis. The required derivatives are described in terms of Lagrange multipliers associated with the value functions.

In this paper, focusing our attention on the dual of a Multistage Stochastic Linear Program (MSLP), we are able to
provide insights into  three important problems for MSLPs: sensitivity analysis,
computation of a sequence of deterministic upper bounds on the optimal value which converges to the optimal value, and use of duality
to solve Dynamic Programming (DP) equations on the dual which are simpler to solve (in the sense that they have convex
cost-to-go functions) than primal DP equations for problems with interstage dependent cost coefficients. Our main contributions are summarized below.

\par {\textbf{Sensitivity analysis of MSLPs.}} We explain how to compute derivatives
of the optimal value, seen as a function of the problem parameters, of a
MSLP in terms of dual optimal solutions. Therefore, the construction of the dual
problem is essential for our approach, contrary to \cite{bonnan}.
With respect to the sensitivity analysis \cite{woza}, in our approach, we do not use value functions directly, which are not known and can only be approximated, but rather construct the dual problem which is solved by an SDDP type algorithm, called
Dual SDDP.

\par {\textbf{Writing Dynamic Programming equations for the dual problem.}} A simple but crucial
ingredient for our developments and subsequent analysis of solution methods for the dual problem of a MSLP
is to write DP equations for that dual problem. We are not aware of another paper with these equations.
However, a similar study was done in \cite{cermics2018}. More precisely, for a stochastic linear
control problem with uncertainty in the right-hand-side, in \cite{cermics2018}, DP equations are
written for the conjugate of the cost-to-go functions and using an SDDP type method for these DP equations,
a sequence of upper bounds on the MSP optimal value is constructed which is the sequence of conjugate of the approximate
first stage cost-to-go functions evaluated at the initial state $x_0$.
Our approach has the advantage of being much simpler: contrary to derivations in \cite{cermics2018} which
require some algebra, our DP equations can be immediately obtained from the dual problem formulation,
this latter being known (given in \cite{shadenrbook} for instance).
On top of that, we relax two assumptions made in \cite{cermics2018}: (a) the relatively complete recourse
assumption of the dual and (b) randomness in the right-hand-side of the constraints only and interstage
independent. The next three paragraphs describe how the scope of (a) and (b) was extended in our
analysis.

\par {\textbf{Dual SDDP for dual problems without relatively complete recourse.}}
In \cite{cermics2018}, it is assumed that the dual problem of the considered MSLP satisfies
an assumption (Assumption 3) stronger than relatively complete recourse. This assumption may not be easy
to check or may not be satisfied (for instance it is not satisfied for the inventory and hydro-thermal
problems considered in Section \ref{numer}). Therefore, it is desirable to extend the scope of Dual
SDDP in such a way that it can still compute a deterministic converging sequence of upper bounds without
this assumption. We present two variants of Dual SDDP that can do that: Dual SDDP with penalizations
and Dual SDDP with feasibility cuts.

\par {\textbf{Dual SDDP for dual problems with all problem data random.}} Our DP equations are
written for problems with uncertainty in all parameters. We explain how to apply Dual SDDP
for such problems that do not satisfy (b) above.

\par {\textbf{Dual SDDP for problems with interstage dependent cost coefficients.}} Finally, we also
relax assumption (b) considering problems having interstage dependent cost coefficients.
Writing DP equations for the corresponding dual problem, we can apply Dual SDDP algorithm to solve
these equations, which, interestingly, have concave cost-to-go functions whereas primal cost-to-go
functions are not convex. This is in sharp contrast with the solution methods proposed so far such as
\cite{bonnan, shaplondorf} which apply SDDP on the primal cost-to-go functions using a Markov
chain approximation of the cost coefficients process.\\

The outline of the paper is the following.
Our building blocks are elaborated in
Section \ref{sec-dual} where we write DP equations for the dual, we explain how to
build upper bounding functions for the cost-to-go functions of the dual using penalizations, and study the dynamics of Lagrange multipliers. Sensitivity analysis of MSLPs is conducted in Section  \ref{sec-sens}
while Dual SDDP and its variants are studied in Section \ref{sec:dualsddp}.
Finally, the results of numerical simulations testing the tools developed on an inventory and an
hydro-thermal problem are presented in Section \ref{numer}. The interested reader can find
and test the code of all implementations and of Primal and Dual SDDP for MSLPs
at {\url{https://github.com/vguigues/Dual_SDDP_Library_Matlab}} and {\url{https://github.com/vguigues/Primal_SDDP_Library_Matlab}}. Proofs are collected in the Appendix.

\mathcal{S}etcounter{equation}{0}
\mathcal{S}ection{Duality of multistage linear stochastic  programs}\label{sec-dual}

\mathcal{S}ubsection{Writing Dynamic Programming equations for the dual}\label{wrdpequ}

Consider the  multistage linear stochastic program
\begin{equation}gin{equation}\label{eq-1}
  \begin{equation}gin{array}{cll}
    \min\limits
    _{x_t\ge 0} & \mbox{\tiny\boldmath$b$}e\left[  \mathcal{S}um_{t=1}^T c_t^\top x_t\right]\\
    {\rm s.t.}& A_1 x_1=b_1,\\
    & B_t x_{t-1}+A_t x_t=b_t,\;t=2,...,T.
   \varepsilonnd{array}
\varepsilonnd{equation}
Here     vectors $c_t=c_t({ x}i_t)\in \mbox{\tiny\boldmath$b$}r^{n_t}$, $b_t =b_t({ x}i_t)\in \mbox{\tiny\boldmath$b$}r^{m_t}$  and matrices $B_t=B_t({ x}i_t)$, $A_t=A_t({ x}i_t)$ are functions of  random process ${ x}i_t\in \mbox{\tiny\boldmath$b$}r^{d_t}$, $t=1,...,T$   (with ${ x}i_1$ being deterministic). We denote by ${ x}i_{[t]}=({ x}i_1,...,{ x}i_t)$ the history of the data process up to time $t$ and by $\mbox{\tiny\boldmath$b$}e_{|{ x}i_{[t]}}$ the corresponding conditional expectation.

The optimization in \varepsilonqref{eq-1} is performed over functions (policies)  $x_t=x_t({ x}i_{[t]})$, $t=1,...,T,$  of the data process  satisfying the feasibility constraints.

The Lagrangian of problem \varepsilonqref{eq-1} is
\begin{equation}gin{equation}\label{lagr}
\begin{equation}gin{array}{l}
 L(x,\pi)=\mbox{\tiny\boldmath$b$}e\left[  \mathcal{S}um_{t=1}^T c_t^\top x_t+\pi_t^\top (b_t-B_t x_{t-1}-A_t x_t)\right]
 \varepsilonnd{array}
\varepsilonnd{equation}
in variables\footnote{Note that since ${ x}i_1$ is deterministic, the first stage decision  $x_1$  is  also deterministic; we write it as $x_1({ x}i_{[1]})$ for uniformity of notation, and similarly for $\pi_1$.}
$x=(x_1({ x}i_{[1]}),\ldots,x_T({ x}i_{[T]}))$ and
$\pi=(\pi_1({ x}i_{[1]}),\ldots,\pi_{T}({ x}i_{[T]}))$
with the convention that $x_0=0$.
Dualization of the feasibility constraints leads to the following dual of problem \varepsilonqref{eq-1} (cf., \cite[Section 3.2.3]{shadenrbook}):
\begin{equation}gin{equation}\label{eq-2}
\begin{equation}gin{array}{cll}
  \max\limits_{\pi} & \mbox{\tiny\boldmath$b$}e\big[  \mathcal{S}um_{t=1}^T     b_t^\top \pi_t\big]\\
    {\rm s.t.}  &A^\top_T\pi_T\le c_T,\\
  &  A_{t -1}^\top \pi_{t-1}+
    \mbox{\tiny\boldmath$b$}e_{|{ x}i_{[t-1]}}\left[ B_{t}^\top \pi_{t}\right ]\le c_{t-1},\;t=2,...,T.
   \varepsilonnd{array}
\varepsilonnd{equation}
The optimization in \varepsilonqref{eq-2} is over policies $\pi_t=\pi_t({ x}i_{[t]})$, $t=1,...,T$.

Unless stated otherwise, we make the  following assumption throughout the paper.
\begin{equation}gin{itemize}
   \item[(A1)]
    The process ${ x}i_1,...,{ x}i_T$ is stagewise independent  (i.e., random vector  ${ x}i_{t+1}$ is independent of ${ x}i_{[t]}$, $t=1,...,T-1$),
    and distribution of ${ x}i_t$ has a finite support, $\{{ x}i_{t 1},\ldots,{ x}i_{t N_t}\}$ with respective probabilities $p_{tj}$,
    $j=1,...,N_t$,  $t=2,...,T$.
 We denote by $A_{t j},B_{t j},c_{t j}, b_{t j}$ the respective  scenarios corresponding to ${ x}i_{t j}$.
 \varepsilonnd{itemize}

 Since  the random process ${ x}i_t$, $t=1,...,T$, has a {\varepsilonm finite} number of realizations (scenarios),  problem \varepsilonqref{eq-1} can be viewed as a large linear program  and \varepsilonqref{eq-2} as its dual. By the standard theory of linear programming we have the following.

 \begin{equation}gin{proposition}
 \label{pr-duality}
 Suppose   that problem \varepsilonqref{eq-1} has a finite optimal value. Then   the optimal values of problems \varepsilonqref{eq-1} and \varepsilonqref{eq-2} are equal to each other and both problems have optimal solutions.
 \varepsilonnd{proposition}

We can write the following dynamic programming equations for the dual problem \varepsilonqref{eq-2}.  At the last stage $t=T$,
given $\pi_{T-1}$ and
${ x}i_{[T-1]}$, we need to solve the following  problem with respect to $\pi_T$:
\begin{equation}gin{equation}\label{eq-3-1}
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_T} &    \mbox{\tiny\boldmath$b$}e[  b_T^\top \pi_T]\\
    {\rm s.t.}&
  A_{T}^\top \pi_{T}
    \le c_{T},\\
    &A_{T -1}^\top \pi_{T-1}+
    \mbox{\tiny\boldmath$b$}e \left[ B_{T}^\top \pi_{T}\right ]\le c_{T-1}.
   \varepsilonnd{array}
\varepsilonnd{equation}
Since ${ x}i_T$  is independent of ${ x}i_{[T-1]}$,  the expectation in \varepsilonqref{eq-3-1} is unconditional with respect to the distribution of ${ x}i_T$.
In terms of scenarios the above problem can be written as
\begin{equation}gin{equation}\label{eq-3a-1}
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_{T 1},\ldots,\pi_{T N_T}} &  \mathcal{S}um\limits_{j=1}^{N_T}   p_{T j}  b_{T j}^\top \pi_{T j}\\
    {\rm s.t.}&
  A_{T j}^\top \pi_{T j}  \le c_{T j},\;j=1,...,N_T,\\
    &A_{T -1}^\top \pi_{T-1}+
  \mathcal{S}um\limits_{j=1}^{N_T}   p_{T j}  B_{T j}^\top
     \pi_{T j} \le c_{T-1}.
   \varepsilonnd{array}
\varepsilonnd{equation}

The optimal value $V_T(\pi_{T-1},{ x}i_{T-1})$ and an optimal   solution\footnote{Note that problem \varepsilonqref{eq-3a-1} may have more than one optimal solution. In case of finite number of scenarios the considered linear program always has a solution provided its optimal value is finite.}
$(\mbox{\boldmath$a$}r{\pi}_{T 1},\ldots,\mbox{\boldmath$a$}r{\pi}_{T N_T})$  of   problem \varepsilonqref{eq-3a-1}  are  functions of vectors  $\pi_{T-1}$ and $c_{T-1}$ and matrix  $A_{T-1}$.
And  so on going backward in time, using  the stagewise independence assumption,  we can write the respective dynamic programming equations for
  $t=T-1,...,2$, as
 \begin{equation}gin{equation}\label{eq-6-1}
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_{t 1},\ldots,\pi_{t N_t}} &     \mathcal{S}um\limits_{j=1}^{N_{t}}  p_{t j}  \left[ b_{t j}^\top \pi_{t j}+ V_{t+1}(\pi_{t j},{ x}i_{t j})\right]\\
    {\rm s.t.}&
  A_{t-1}^\top \pi_{t-1}+
  \mathcal{S}um\limits_{j=1}^{N_{t}} p_{tj} B_{t j}^\top  \pi_{t j}  \le c_{t-1},
   \varepsilonnd{array}
\varepsilonnd{equation}
with $V_t(\pi_{t-1},{ x}i_{t-1})$ being the optimal value  of problem \varepsilonqref{eq-6-1}.
Finally at the first stage the following problem should be solved
 \begin{equation}gin{equation}\label{eq-7-1}
  \max_{\pi_{1}}      b_1^\top \pi_1+ V_{2}(\pi_1,{ x}i_1) .
   \varepsilonnd{equation}

These dynamic  programming equations can be compared with the
dynamic  programming equations for primal problem \varepsilonqref{eq-1}, where the respective cost-to-go (value) function  $Q_t(x_{t-1},{ x}i_{t j})$, $j=1,...,N_t$,  is  given by the optimal value of
\begin{equation}gin{equation}\label{primal}
\begin{equation}gin{array}{cll}
\min\limits_{x_t\ge 0} & c_{t j}^\top x_t+{\cal Q}_{t+1}(x_t)\\
{\rm s.t.} & B_{t j} x_{t-1}+A_{t j} x_t=b_{t j},
\varepsilonnd{array}
\varepsilonnd{equation}
with
\[
{\cal Q}_{t+1}(x_t):=\mbox{\tiny\boldmath$b$}e[ Q_{t+1}(x_t,{ x}i_{t+1})]=\mathcal{S}um_{j=1}^{N_{t}} p_{t+1 j} Q_{t+1}(x_t,{ x}i_{t+1 j}).
\]

Let us make the following observations  about   the dual problem.
\begin{equation}gin{itemize}
  \item [(i)]
  Unlike in  the primal problem,
  the optimization (maximization) problems \varepsilonqref{eq-3a-1} and \varepsilonqref{eq-6-1}   do not decompose into separate problems with respect to each $\pi_{t j}$ and should be solved as one linear program with respect to $(\pi_{t 1},...,\pi_{t N_t})$.
  \item [(ii)]
 The value function   $V_t(\pi_{t-1},{ x}i_{t-1})$   is a concave function of $\pi_{t-1}$.
 \item [(iii)]
If $A_t$ and $c_t$, $t=2,...,T$,  are deterministic, then $V_t(\pi_{t-1})$ is only a function  of $\pi_{t-1}$.
\varepsilonnd{itemize}

\mathcal{S}ubsection{Relatively complete recourse}\label{sec:rcr}

The following definition of Relatively Complete Recourse (RCR) is applied to the dual problem.    Recall that we assume that the set of possible realizations (scenarios) of the data process is finite.
\begin{equation}gin{definition}
\label{def-rcr}
We say that a sequence $\mbox{\boldmath$a$}r{\pi}_t$, $t=1,...,T$, is generated by the forward (dual)  process if $\mbox{\boldmath$a$}r{\pi}_1\in \mbox{\tiny\boldmath$b$}r^{m_1}$ and   for    $\pi_{t-1}=\mbox{\boldmath$a$}r{\pi}_{t-1}$,   $t=2,...,T$, going forward in time,  $\mbox{\boldmath$a$}r{\pi}_t$  coincides with some $\pi_{t j}$, $j=1,...,N_t$, where
$\pi_{t 1},\ldots,\pi_{t N_t}$ is a feasible solution
of the respective dynamic program - program  \varepsilonqref{eq-6-1} for  $t=2,...,T-1$,  and program \varepsilonqref{eq-3a-1} for $t=T$.
We say that the dual problem \varepsilonqref{eq-2} has {\rm Relatively
Complete Recourse (RCR)} if    at every stage $t=2,...,T$,  for  any  generated $\pi_{t-1}$  by the forward process,    the respective dynamic program has a  feasible  solution at stage $t$ for every  realization of the random data.
\varepsilonnd{definition}

Without  RCR   it could  happen that   $V_t(\pi_{t-1},{ x}i_{t-1})=-\infty$ for a generated  $\pi_{t-1}$   and ${ x}i_{t-1}={ x}i_{t-1 j}$. Unfortunately, it could happen that  the dual problem does not have the RCR property even if the primal problem has it. This could happen even in the two stage case.
One way to deal with     the  problem of absence of RCR  in numerical procedures is to use feasibility cuts, we will discuss this later. Another way  is the following penalty approach which will be used in Section \ref{sec:dualsddp}.

The infeasibility of problem \varepsilonqref{eq-3a-1} can happen because of its last constraint.  In order to deal with this, consider the following relaxation of problem \varepsilonqref{eq-3a-1}:
\begin{equation}gin{equation}\label{eq-r1}
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_{T 1},...,\pi_{T N_T},\!\!,\,{ z}eta_T\ge 0} &    \mathcal{S}um\limits_{j=1}^{N_T}   p_{T j}  b_{T j}^\top \pi_{T j}-v_T^\top { z}eta_T\\
    {\rm s.t.}&
  A_{T j}^\top \pi_{T j}  \le c_{T j},\;j=1,...,N_T,\\
    &A_{T -1}^\top \pi_{T-1}+
     \mathcal{S}um\limits_{j=1}^{N_T}   p_{T j}  B_{T j}^\top \pi_{T j} \le c_{T-1}+{ z}eta_T,
   \varepsilonnd{array}
\varepsilonnd{equation}
where $v_T$ is a  vector with positive components.
We have that problem \varepsilonqref{eq-r1} is always feasible
and hence its  optimal value $\tilde{V}_T(\pi_{T-1},{ x}i_{T-1})>-\infty$.
We also have that
\begin{equation}gin{equation}\label{valineq}
\tilde{V}_T(\pi_{T-1},{ x}i_{T-1})\ge V_T(\pi_{T-1},{ x}i_{T-1}),
\varepsilonnd{equation}
 with   the equality
 holding  if ${ z}eta_T=0$  in the optimal solution of \varepsilonqref{eq-r1}.  If
$V_T(\pi_{T-1}$,${ x}i_{T-1})$ is finite, this equality holds if the components of vector $v_T$ are large enough.

Similarly, problems \varepsilonqref{eq-6-1} can be relaxed to
\begin{equation}gin{equation}\label{eq-r2}
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_{t 1},\ldots,\pi_{t N_t}, { z}eta_t \geq 0} &  \mathcal{S}um\limits_{j=1}^{N_{t}}  p_{t j}  \left[ b_{t j}^\top
  \pi_{t j}+ \tilde{V}_{t+1}(\pi_{t j},{ x}i_{t j})\right]-v_t^\top { z}eta_t\\
    {\rm s.t.}&
  A_{t-1}^\top \pi_{t-1}+
  \mathcal{S}um\limits_{j=1}^{N_{t}} p_{tj} B_{t j}^\top  \pi_{t j}  \le c_{t-1}+{ z}eta_t,
   \varepsilonnd{array}
\varepsilonnd{equation}
with vector  $v_t$ having positive components.
In that way, the infeasibility problem is avoided and the  obtained  value gives an  upper bound for the optimal value of the dual problem.   Note that for sufficiently large vectors $v_t$ this upper bound coincides with the optimal value of the dual problem.

\mathcal{S}ubsection{Dynamics of Lagrange multipliers}
\label{sec-dyn}

Let us consider for the moment the two stage setting, i.e., $T=2$. The primal problem can be written as

\begin{equation}gin{equation}\label{lag-1}
    \min\limits_{x_1\ge 0}   c_1^\top x_1+ \mbox{\tiny\boldmath$b$}e\left[ Q(x_1,{ x}i_2)\right]\;
    {\rm s.t.}\; A_1 x_1=b_1,
 \varepsilonnd{equation}
where $Q(x_1,{ x}i_2)$ is the optimal value of the second stage problem
\begin{equation}gin{equation}\label{lag-2}
    \min\limits_{x_2\ge 0}      c_2({ x}i_2)^\top x_2\;
    {\rm s.t.}\;  B_2({ x}i_2)  x_{1}+A_2({ x}i_2)  x_2 =b_2({ x}i_2).
\varepsilonnd{equation}
The Lagrangian of problem \varepsilonqref{lag-2} is
\[
L(x_1,x_2,\lambda,{ x}i_2)=c_2({ x}i_2)^\top x_2+\lambda^\top (b_2({ x}i_2)
- B_2({ x}i_2)  x_{1}-A_2({ x}i_2)  x_2).
\]
In the dual form,   $Q(x_1,{ x}i_{2 j})$ is given by the optimal value of the problem

 \begin{equation}gin{equation}\label{lag-3}
 \max\limits_{\lambda_j }  \; (b_{2 j}-B_{2 j} x_1)^\top \lambda_j   \;
  {\rm s.t.}
 \;  c_{2 j} -    A_{2 j}^\top \lambda_j   \ge 0.
\varepsilonnd{equation}
We have that if $x_1=\mbox{\boldmath$a$}r{x}_1$ is an optimal solution of the first stage problem, then optimal Lagrange multipliers $\pi_{2 j}$ are given by the optimal solution of problem \varepsilonqref{lag-3}.

This can be extended to the multistage setting of problem \varepsilonqref{eq-1} (recall that the stagewise independence condition is assumed). At the last stage $t=T$, given optimal solution $\mbox{\boldmath$a$}r{x}_{T-1}$, the following problem should be solved
\begin{equation}gin{equation}\label{lag-4}
    \min\limits_{x_T\ge 0}     c_T({ x}i_T)^\top x_T \;\;
    {\rm s.t.}\;\;  B_T({ x}i_T) \mbox{\boldmath$a$}r{x}_{T-1}+A_T({ x}i_T) x_T=b_T({ x}i_T).
\varepsilonnd{equation}
For a realization  ${ x}i_T={ x}i_{T j}$,  the dual of problem \varepsilonqref{lag-4} is the problem
 \begin{equation}gin{equation}\label{lag-5}
 \max\limits_{\lambda_j } \;  (b_{T j}-B_{T j} \mbox{\boldmath$a$}r{x}_{T-1})^\top \lambda_j   \;\; {\rm s.t.}  \;\;  c_{T j} -  A_{T j}^\top \lambda_j   \ge 0.
\varepsilonnd{equation}
We then have that $\pi_{T j}$ are given by the optimal solution of problem \varepsilonqref{lag-5}.

At stage $t=T-1$, given optimal solution $\mbox{\boldmath$a$}r{x}_{T-2}$,  the following problem is supposed to  be solved (see  \varepsilonqref{primal})
\begin{equation}gin{equation}\label{lag-6}
\begin{equation}gin{array}{cll}
    \min\limits_{x_{T-1}\ge 0}  &   c_{T-1}({ x}i_{T-1})^\top x_{T-1}+{\cal Q}_T(x_{T-1})\\
    {\rm s.t.} &  A_{T-1}({ x}i_{T-1}) x_{T-1}=b_{T-1}({ x}i_{T-1})-B_{T-1}({ x}i_{T-1}) \mbox{\boldmath$a$}r{x}_{T-2}.
    \varepsilonnd{array}
\varepsilonnd{equation}
   We have that ${\cal Q}_T(\cdot)$  is a convex piecewise linear function. Therefore for every realization ${ x}i_{T-1}={ x}i_{T-1 j}$  it is possible to represent \varepsilonqref{lag-6} as a linear program and hence to write its dual.  The optimal  Lagrange multipliers of that dual give the corresponding
Lagrange multipliers $\pi_{T-1 j}$.  And so on for other  stages going backward in time.  That is, we have the following.
\begin{equation}gin{remark}\label{remarkcomputedualmpsddp}  If $(\mbox{\boldmath$a$}r{x}_1,...,\mbox{\boldmath$a$}r{x}_T({ x}i_{[T]}))$ is an optimal solution of the primal problem, then for $x_{t-1}=\mbox{\boldmath$a$}r{x}_{t-1}$ the Lagrange multiplier $\pi_{t j}$ is given by the respective  Lagrange multiplier of problem \varepsilonqref{primal}.
\varepsilonnd{remark}

We also refer to \cite{higlesenaor, rockdualsol} for the dynamics
of dual solutions to stochastic programs.

\mathcal{S}etcounter{equation}{0}
\mathcal{S}ection{Sensitivity analysis}
\label{sec-sens}
In this section we discuss an application of the duality analysis to a study of sensitivity of the optimal value to small perturbations of the involved parameters.

\mathcal{S}ubsection{General case}
Suppose now that the data $c_t({ x}i_t,\theta),  b_t({ x}i_t,\theta), B_t({ x}i_t,\theta)$, $A_t({ x}i_t$, $\theta)$ of problem \varepsilonqref{eq-1} also depend  on parameter vector $\theta\in \mbox{\tiny\boldmath$b$}r^k$. Denote by
$\vartheta(\theta)$  the optimal value of the parameterized problem   \varepsilonqref{eq-1}  considered as a function of $\theta$, and by ${\mathfrak S}(\theta)$ and ${\mathfrak D}(\theta)$ the sets of optimal solutions of the respective primal and dual problems. Recall that the sets ${\mathfrak S}(\theta)$ and ${\mathfrak D}(\theta)$ are nonempty provided the optimal value $\vartheta(\theta)$ is finite.
Let $L(x,\pi,\theta)$ be the corresponding Lagrangian (see \varepsilonqref{lagr}) considered  as a function of $\theta$.
Then
we have the following formula for the directional derivatives of the optimal value function (e.g.,  \cite[Proposition 4.27]{Bonnans2000}).

\begin{equation}gin{proposition}
\label{pr-valsens}
Suppose   that    the data functions are continuously differentiable functions of $\theta$, and   for a given $\theta=\mbox{\boldmath$a$}r{\theta}$ the optimal value $\vartheta(\mbox{\boldmath$a$}r{\theta})$ is finite and the sets ${\mathfrak S}(\mbox{\boldmath$a$}r{\theta})$ and ${\mathfrak D}(\mbox{\boldmath$a$}r{\theta})$ of optimal solutions   are bounded. Then
\begin{equation}gin{equation}\label{sen-1}
 \vartheta'(\mbox{\boldmath$a$}r{\theta},h)=\max_{\pi\in {\mathfrak D}(\mbox{\boldmath$a$}r{\theta})}\min_{x\in {\mathfrak S}(\mbox{\boldmath$a$}r{\theta})} h^\top \nabla_\theta L(x,\pi,\mbox{\boldmath$a$}r{\theta}).
\varepsilonnd{equation}
In particular if ${\mathfrak S}(\mbox{\boldmath$a$}r{\theta})=\{\mbox{\boldmath$a$}r{x}\}$ and ${\mathfrak D}(\mbox{\boldmath$a$}r{\theta})=\{\mbox{\boldmath$a$}r{\pi}\}$
are singletons, then $\vartheta(\cdot)$  is differentiable at $\mbox{\boldmath$a$}r{\theta}$ and
\begin{equation}gin{equation}\label{sen-2}
 \nabla \vartheta(\mbox{\boldmath$a$}r{\theta})= \nabla_\theta L(\mbox{\boldmath$a$}r{x},\mbox{\boldmath$a$}r{\pi},\mbox{\boldmath$a$}r{\theta}).
\varepsilonnd{equation}
\varepsilonnd{proposition}

Next, as an example, we consider the sensitivity analysis of an inventory model.

\mathcal{S}ubsection{Application to an inventory model}
\label{sec-invent}

Consider the inventory model
\begin{equation}gin{equation}\label{inv-1}
\begin{equation}gin{array}{cll}
 \min & \mathbb{E}\left[\displaystyle \mathcal{S}um_{t=1}^T
 a_t( y_t - x_{t-1} ) + g_t (\mathcal{D}_t - y_t)_{+} + h_t (y_t - \mathcal{D}_t)_{+} \right]\\
 {\rm s.t.}& x_{t} = y_t -\mathcal{D}_t, y_t \geq x_{t-1},t=1\ldots,T.
 \varepsilonnd{array}
\varepsilonnd{equation}
Here  $\mathcal{D}_1,...,\mathcal{D}_T$ is a (random) demand process,   $a_t,g_t,h_t$ are the ordering, back-order penalty and holding costs per unit, respectively,
$x_t$ is the inventory level and $y_t-x_{t-1}$ is the order quantity at time $t$,
 the initial inventory level $x_0$ is given. We refer to   \cite{Zipk} for a thorough discussion of that model.
Note that $\mathcal{D}_t$
is a random variable whereas
$d_t$ stands for a particular realization.
We assume that $g_t>a_t\ge 0$, $h_t> 0$, $t=1,...,T$.

In the classical setting the demand process is assumed to be stagewise independent, i.e., $\mathcal{D}_{t+1}$ is assumed to be independent of $\mathcal{D}_{[t]}=(\mathcal{D}_1,...,\mathcal{D}_t)$ for $t=1,...,T-1$.

In order to capture the autocorrelation structure of the demand process it is tempting  to model it as, say first order, autoregressive process
$
 \mathcal{D}_t=\mu+\phi \mathcal{D}_{t-1}  +{\varepsilonpsilon}_t,
$
where errors  ${\varepsilonpsilon}_t$ are assumed to be a sequence  i.i.d (independent identically distributed)  random variables. However this approach may result in some of the realizations of the demand process to be negative, which of course does not make sense. One way to deal with this is to make the transformation  $Y_t:=\log \mathcal{D}_t$ and to model $Y_t$  as an autoregressive process. A problem with this approach is  that it  leads to nonlinear equations for the original process $\mathcal{D}_t$,  which makes it difficult to use in the numerical algorithms discussed below.

We assume that
the demand  is modeled as the following multiplicative autoregressive process
\begin{equation}gin{equation}
\label{inv-2-1bis}
\mathcal{D}_t = {\varepsilonpsilon}_{t}( \phi \mathcal{D}_{t-1} + \mu), \ t = 1,...,T,
\varepsilonnd{equation}
where $\phi\in (0,1)$, $\mu\ge 0$ are   parameters and $\mathcal{D}_0\ge 0$ is given. The errors ${\varepsilonpsilon}_t$ are i.i.d  with  log-normal distributions having means and standard deviations given by $\mbox{\tiny\boldmath$b$}e[{\varepsilonpsilon}_t]=1$
and $\text{Var}({\varepsilonpsilon}_t)=\mathcal{S}igma^2>1$, respectively.
This guarantees that all realizations of the demand process are positive.   It is possible to view \varepsilonqref{inv-2-1bis} as a linearization of the log-transformed process
$\log \mathcal{D}_t$  (cf., \cite{sha2013}).
See Section \ref{statpr} for a discussion of statistical properties of the  process \varepsilonqref{inv-2-1bis}.

The process \varepsilonqref{inv-2-1bis} involves parameters $\phi$ and $\mu$ which are supposed to be estimated from the data. As such, these parameters are subject to estimation errors. This raises the question of sensitivity of the optimal value $\vartheta=\vartheta(\phi,\mu)$  of the corresponding problem \varepsilonqref{inv-1} viewed as a function of $\phi$ and $\mu$. To this end, we investigate the calculation of the derivatives
$\partial \vartheta(\phi,\mu)/\partial \phi$ and $\partial \vartheta(\phi,\mu)/\partial \mu$.  With these derivatives at hand, asymptotic distributions of the estimates of $\phi$ and $\mu$ can be translated into the asymptotics of the optimal value in a straightforward way by application of the Delta Theorem.  We refer to Section \ref{invnumer} for the corresponding numerical experiments.

\mathcal{S}ubsubsection{Properties of
 the multiplicative autoregressive process}
 \label{statpr}

Consider the multiplicative autoregressive process \varepsilonqref{inv-2-1bis}.
Note that under the specified conditions the demand process is not stationary. Indeed,
since the errors ${\varepsilonpsilon}_t$ are i.i.d and $\mbox{\tiny\boldmath$b$}e[{\varepsilonpsilon}_t]=1$  we have  that
$
\mbox{\tiny\boldmath$b$}e[\mathcal{D}_t ]=   \phi \mbox{\tiny\boldmath$b$}e[\mathcal{D}_{t-1}] + \mu
$
and
\begin{equation}gin{equation}
\label{inv-3}
\begin{equation}gin{array}{lll}
{\rm Var}(\mathcal{D}_t)&=&\mbox{\tiny\boldmath$b$}e\left[{\rm Var}\big ( \varepsilon_{t}( \phi \mathcal{D}_{t-1} + \mu)|\mathcal{D}_{t-1}\big )\right]+  {\rm Var}  \left[ \mbox{\tiny\boldmath$b$}e (\varepsilon_{t}( \phi \mathcal{D}_{t-1} + \mu)|\mathcal{D}_{t-1})\right] \\
&=&\mbox{\tiny\boldmath$b$}e\left[ \mathcal{S}igma^2 ( \phi \mathcal{D}_{t-1} + \mu)^2\right]+{\rm Var}( \phi \mathcal{D}_{t-1} + \mu)\\
&=&\mathcal{S}igma^2 \mbox{\tiny\boldmath$b$}e\left[( \phi \mathcal{D}_{t-1} + \mu)^2\right]+\phi^2 {\rm Var}(\mathcal{D}_{t-1} ).
\varepsilonnd{array}
\varepsilonnd{equation}
It follows that $\mbox{\tiny\boldmath$b$}e[\mathcal{D}_t]$ converges to $\mu/(1-\phi)$ as $t\to \infty$.
 Suppose, for example,   that $\mu=0$. Then
 $
\mathcal{D}_t = \varepsilon_{t}  \phi \mathcal{D}_{t-1}=\mathcal{D}_0\phi^t\prod_{\tau=1}^t \varepsilon_\tau, \ t = 1,...,T,
$ $\mbox{\tiny\boldmath$b$}e[\mathcal{D}_t]=\mathcal{D}_0 \phi^t\to 0$, and
 $
{\rm Var}(\mathcal{D}_t )=\mathcal{D}_0^2\phi^{2t}[(1+\mathcal{S}igma^2)^t-1].
$
Therefore if $ \phi^2(1+\mathcal{S}igma^2)<1$, then
${\rm Var}(\mathcal{D}_t )\to 0$; and  if $ \phi^2(1+\mathcal{S}igma^2)> 1$, then
${\rm Var}(\mathcal{D}_t )\to \infty$ provided $\mathcal{D}_0>0$.

\mathcal{S}ection{Dual SDDP}\label{sec:dualsddp}

In this section, using the results of Section \ref{sec-dual}, we discuss  an adaptation of the cutting planes approach  for the approximation  of the value functions of the dual problem, similar to the standard SDDP method
and called Dual SDDP. The interested reader
can find the implementation
of  Primal SDDP and all variants of Dual SDDP described in this section at
{\url{https://github.com/vguigues/Dual_SDDP_Library_Matlab}} and {\url{https://github.com/vguigues/Primal_SDDP_Library_Matlab}}.

We will make the following assumption.
\begin{equation}gin{itemize}
\item[(A2)] Primal problem \varepsilonqref{eq-1}
satisfies the RCR assumption.
\varepsilonnd{itemize}

We first consider the case where only $b_t$
and $B_t$ are random in ${ x}i_t$.

\mathcal{S}ubsection{Dual SDDP for problems with uncertainty in $b_t$ and $B_t$}

In Dual SDDP, concave value functions $V_t,t=2,\ldots,T$, are approximated at the end of  iteration $k$ by polyhedral upper bounding functions
$V_t^k$ given by:
\begin{equation}gin{equation}\label{cutvtkdsddp}
V_t^k( \pi_{t-1}) = \displaystyle \min_{0 \leq i \leq k}
{\overline{\theta}}_t^i + \langlegle {\overline{\begin{equation}ta}}_t^i , \pi_{t-1} {\ranglegle}ngle
\varepsilonnd{equation}
where
${\overline{\theta}}_t^i$, ${\overline{\begin{equation}ta}}_t^i$
are coefficients whose computation is detailed below. The algorithm uses valid upper bounds on the norm of dual optimal solutions:

\begin{equation}gin{lemma}\label{boundedpi}
Suppose that the optimal value of primal problem \varepsilonqref{eq-1} is finite  and that there is ${\hat x}>0$
feasible for primal problem
\varepsilonqref{eq-1}.
Then for every $t=1,\ldots,T$,
we can find
$\underline{\pi}_t, \overline{\pi}_t
\in \mathbb{R}^{m_t}$
such that dual problem
\varepsilonqref{eq-6-1}  is unchanged (i.e., has the
same optimal value) adding box
constraints
$\underline{\pi}_t \leq \pi_t \leq  \overline{\pi}_t$.

\varepsilonnd{lemma}

Recall that it is assumed that the number of scenarios is finite and hence problem \varepsilonqref{eq-1} can be viewed as a large linear program. The assumption of existence of feasible ${\hat x}>0$ means that   problem  \varepsilonqref{eq-1} possesses a feasible solution with all components being strictly positive.
If moreover  the equality constraints of problem \varepsilonqref{eq-1} are linearly independent, then this strict feasibility condition   implies that the set of optimal solutions of the dual problem (i.e., the set of
Lagrange multipliers) is bounded. On the other hand,  in the above lemma the linear independence condition is not assumed. A proof of Lemma \ref{boundedpi} and a way to obtain the corresponding
bounds $\underline{\pi}_t, \overline{\pi}_t$
can be found in the Appendix.

As mentioned earlier, a difficulty to solve the dual
problem with an SDDP type method is that
RCR may not be satisfied by the dual problem, even if RCR holds for the primal.
We propose two variants of Dual SDDP
to solve the Dual problem even if RCR does not hold for the dual:
Dual SDDP with penalizations and Dual SDDP with feasibility cuts.

\par {\textbf{Dual SDDP with penalizations.}} Dual SDDP with penalizations is based on the developments of
Section \ref{sec:rcr}. It introduces slack variables in the constraints which may become infeasible
for some past decisions in the subproblems solved in the forward passes of Dual SDDP.
Slack variables are penalized in the objective function with
sequences $(v_{t k})_k$  of positive penalizing coefficients.
Therefore, all subproblems solved in forward and backward passes of this variant of Dual SDDP, called
Dual SDDP with penalizations, are always feasible
and at iteration $k$, the method builds
polyhedral upper bounding function  $V_t^k$ for $V_t$ of form \varepsilonqref{cutvtkdsddp} (see Proposition \ref{proofubfeassddp}). Similarly to SDDP, trial points are generated in a forward pass and cuts for $V_t$
are computed in a backward pass. The detailed Dual SDDP method with penalizations
is as follows.\\

\par {\textbf{Initialization.}} For $t=2,\ldots,T,$ take for ${V}_t^0$
an affine upper bounding function for $V_t$
and $V_{T+1}^0 \varepsilonquiv 0$. Set iteration counter $k$ to 1.\\

\par {\textbf{Step 1: forward pass of iteration $k$ (computation of dual trial points).}}
For the first stage of the forward pass, we compute an optimal solution
$\pi_1^k$ of
\begin{equation}gin{equation}\label{firststagepenalized}
V^{k-1} =
\begin{equation}gin{array}{l}
\displaystyle \max_{\pi_1} \; b_1^\top \pi_1 + V_{2}^{k-1}(\pi_1)\\
{\underline{\pi}}_1 \leq \pi_1 \leq {\overline{\pi}}_1.
\varepsilonnd{array}
\varepsilonnd{equation}
Recall that the optimal value of the first stage problem
does not change adding box constraints
${\underline{\pi}}_1 \leq \pi_1 \leq {\overline{\pi}}_1$  for appropriate values ${\underline{\pi}}_1$ and ${\overline{\pi}}_1$.
The introduction of these box constraints ensures that
the optimal value of \varepsilonqref{firststagepenalized} (which is an
{\varepsilonm{approximate}} first stage problem
due to the approximation of
$V_2$ by $V_2^{k-1}$) is finite for all iterations.

\par For stage $t=2,\ldots,T-1$, given $\pi_{t-1}^k$, we compute an optimal solution of
\begin{equation}gin{equation}\label{forwarddualsddppen}
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_{t 1},\ldots,\pi_{t N_t}, { z}eta_t \geq 0} &  \mathcal{S}um\limits_{j=1}^{N_{t}}  p_{t j}  \left[ b_{t j}^\top
  \pi_{t j}+ {V}_{t+1}^{k-1}(\pi_{t j})\right]-v_{t k}^\top { z}eta_t\\
    {\rm s.t.}&
  A_{t-1}^\top \pi_{t-1}^k +
  \mathcal{S}um\limits_{j=1}^{N_{t}} p_{tj} B_{t j}^\top  \pi_{t j}  \le c_{t-1}+{ z}eta_t,\\
  & {\underline{\pi}}_t \leq \pi_{t j} \leq {\overline{\pi}}_t.
   \varepsilonnd{array}
\varepsilonnd{equation}
\par An optimal solution of the problem above has $N_t$ components
$(\pi_{t 1}, \pi_{t 2},\ldots,\pi_{t N_t})$
for $\pi_t$. We generate a realization $\tilde { x}i_t^k$ of
${ x}i_t^k \mathcal{S}im { x}i_t$ independently of previous realizations
$\tilde { x}i_2^1,\ldots$, $\tilde { x}i_{T-1}^1$,$\ldots$,
$\tilde { x}i_2^k,\ldots,\tilde { x}i_{t-1}^k$, and take $\pi_{t}^k = \pi_{t j_t(k)}$
where index $j_t(k)$ satisfies
$\tilde { x}i_{t}^k = { x}i_{t j_t(k)}$.\\

\par {\textbf{Step 2: backward pass of iteration $k$ (computation of new cuts)}}. We first compute a new cut for $V_T$. Let
$(\alpha, \delta, \overline{\Psi}, \underline{\Psi})$ be an optimal solution of\footnote{We suppressed the dependence of the optimal solution on $T$ and $k$ to alleviate notation.}
\begin{equation}gin{equation}\label{defVTk}
\begin{equation}gin{array}{l}
\displaystyle \min_{\alpha, \delta, \overline{\Psi}, {\underline{\Psi}}} \;
\delta^{\top} (c_{T-1} - A_{T-1}^{\top} \pi_{T-1}^k)  + c_T^{\top}  \mathcal{S}um_{j=1}^{N_T} \alpha_j
+ \mathcal{S}um_{j=1}^{N_T} {\overline{\Psi}}_j^\top {\overline \pi}_T
- \mathcal{S}um_{j=1}^{N_T} {\underline{\Psi}}_j^\top {\underline \pi}_T \\
A_T \alpha_j + p_{T j} B_{T j} \delta -
{\underline{\Psi}}_j +{\overline{\Psi}}_j = p_{T j} b_{T j},\;j=1,\ldots,N_T,\\
0 \leq \delta \leq v_{T k}, \alpha_j, {\underline{\Psi}}_j, {\overline{\Psi}}_j \geq 0\;j=1,\ldots,N_T.
\varepsilonnd{array}
\varepsilonnd{equation}
\par The new cut for $V_T$ has  coefficients
given by
$$
\overline \theta_T^k =
\delta^{\top} c_{T-1}+ c_T^{\top}  \mathcal{S}um_{j=1}^{N_T} \alpha_j
+ \mathcal{S}um_{j=1}^{N_T} {\overline{\Psi}}_j^\top {\overline \pi}_T
- \mathcal{S}um_{j=1}^{N_T} {\underline{\Psi}}_j^\top {\underline \pi}_T,\;
\overline \begin{equation}ta_T^k = -A_{T-1}\delta.
$$
\par For $t=T-1,\ldots,2$, compute an optimal solution
$(\delta,\nu,{\overline \Psi}, {\underline \Psi})$ of
\begin{equation}gin{equation}\label{defVtk}
\begin{equation}gin{array}{l}
\displaystyle \min_{\delta,\nu,{\overline \Psi}, {\underline \Psi}} \;
\delta^{\top}{\cal B}ig[
c_{t-1} - A_{t-1}^{\top} \pi_{t-1}^k {\cal B}ig]
+ \mathcal{S}um_{i=0}^{k} {\overline{\theta}}_{t+1}^i \mathcal{S}um_{j=1}^{N_t} \nu_i(j) +
\mathcal{S}um_{j=1}^{N_t} {\overline \Psi}_j^\top {\overline \pi}_t
-  \mathcal{S}um_{j=1}^{N_t} {\underline \Psi}_j^\top {\underline \pi}_t\\
\displaystyle
p_{t j} B_{t j} \delta - \mathcal{S}um_{i=0}^k \nu_i(j)  {\overline{\begin{equation}ta}}_{t+1}^i - {\underline \Psi}_j
+ {\overline \Psi}_j = p_{t j} b_{t j},\;j=1,\ldots,N_t,\\
\displaystyle \mathcal{S}um_{i=0}^k \nu_i( j ) = p_{t j}, {\underline \Psi}_j, {\overline \Psi}_j \geq 0,\;j=1,\ldots,N_t,\\
\nu_0,\ldots,\nu_k \geq 0,0 \leq \delta \leq v_{t k},
\varepsilonnd{array}
\varepsilonnd{equation}
\par and the cut coefficients
$$
\overline{\theta}_t^k = \delta^{\top} c_{t-1}
+ \mathcal{S}um_{i=0}^{k} {\overline{\theta}}_{t+1}^i \mathcal{S}um_{j=1}^{N_t} \nu_i(j) +
\mathcal{S}um_{j=1}^{N_t} {\overline \Psi}_j^\top {\overline \pi}_t
-  \mathcal{S}um_{j=1}^{N_t} {\underline \Psi}_j^\top {\underline \pi}_t,\;
{\overline{\begin{equation}ta}}_t^k = -A_{t-1} \delta.
$$
\par {\textbf{Step 3:}} Do $k \leftarrow k+1$ and go to Step 1.\\

\par The validity of the cuts computed in the backward pass of Dual SDDP with penalizations is shown in Proposition \ref{proofubfeassddp}.
\begin{equation}gin{proposition}\label{proofubfeassddp}
Consider Dual SDDP algorithm with penalizations.
Let Assumptions (A1) and (A2) hold.
Then for every $t=2,\ldots,T$, the sequence $V_t^k$ is a nonincreasing sequence
of upper bounding functions for $V_t$, i.e., for every $k \geq 1$ we have
$V_t \leq {V}_t^k \leq V_t^{k-1}$
 and therefore $(V^k)$ (recall that $V^{k-1}$ is the optimal value of \varepsilonqref{firststagepenalized}) is a nonincreasing deterministic
 sequence of upper bounds on the optimal value of \varepsilonqref{eq-1}.
\varepsilonnd{proposition}
To understand the effect of the sequence of penalizing parameters $(v_{t k})$
on Dual SDDP with penalizations, we define the following Dynamic Programming equations (see also
Lemma \ref{lem-conver} in the Appendix):
\begin{equation}gin{equation}\label{dualdpesimp1bb}
V_T^{\gamma}(\pi_{T-1}) =
\left\{
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_{T 1},...,\pi_{T N_T},\!\!,\,{ z}eta_T\ge 0} &    \mathcal{S}um\limits_{j=1}^{N_T}   p_{T j}  b_{T j}^\top \pi_{T j}-\gamma {\textbf{e}}^\top { z}eta_T\\
    {\rm s.t.}&
  A_{T j}^\top \pi_{T j}  \le c_{T j},\;j=1,...,N_T,\\
    &A_{T -1}^\top \pi_{T-1}+
     \mathcal{S}um\limits_{j=1}^{N_T}   p_{T j}  B_{T j}^\top \pi_{T j} \le c_{T-1}+{ z}eta_T,
   \varepsilonnd{array}
   \right.
   \varepsilonnd{equation}
for $t=2,\ldots,T-1$:
\begin{equation}gin{equation}\label{dualdpesimp2bb}
V_t^{\gamma}(\pi_{t-1}) =
\left\{
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_{t 1},\ldots,\pi_{t N_t}, { z}eta_t \geq 0} &  \mathcal{S}um\limits_{j=1}^{N_{t}}  p_{t j}  \left[ b_{t j}^\top
  \pi_{t j}+ V_{t+1}^{\gamma}(\pi_{t j})\right]-\gamma {\textbf{e}}^\top { z}eta_t\\
    {\rm s.t.}&
  A_{t-1}^\top \pi_{t-1}+
  \mathcal{S}um\limits_{j=1}^{N_{t}} p_{tj} B_{t j}^\top  \pi_{t j}  \le c_{t-1}+{ z}eta_t,
   \varepsilonnd{array}
\right.
\varepsilonnd{equation}
and we define the first stage problem
 \begin{equation}gin{equation}\label{eq-7simpbb}
  \max_{\pi_1} \;\;      \pi_1^\top  b_1 +  V_{2}^{\gamma}(\pi_1),
\varepsilonnd{equation}
where {\textbf{e}} is a vector of ones and $\gamma$ is a positive real number.
As we will see below, $V_t^{\gamma}$ can be seen as an upper bounding concave approximation
of $V_t$ which gets ``closer'' to $V_t$ when $\gamma$ increases.
For inventory problem \varepsilonqref{inv-1}, it is easy to see that
functions $V_t$ in DP equations \varepsilonqref{eq-3a-1}, \varepsilonqref{eq-6-1}, \varepsilonqref{eq-7-1} and functions $V_t^{\gamma}$ in DP equations
\varepsilonqref{dualdpesimp1bb}, \varepsilonqref{dualdpesimp2bb}, \varepsilonqref{eq-7simpbb} (obtained
using in these equations data $c_t, b_t$ $A_t$, $B_t$, corresponding to the inventory problem) are only functions of one-dimensional state
variable $\pi_{t-1}$. Therefore, Dynamic Programming can be used to solve
these Dynamic Programming equations and obtain
good approximations of functions $V_t$ and $V_t^{\gamma}$. To obtain these approximations,
we need to obtain approximations of the domains of functions  $V_t$
and compute approximations of these functions on a set of points in that domain.
To observe the impact of penalizing term $\gamma$ on $V_t^{\gamma}$, we run Dynamic Programming both on DP equations \varepsilonqref{eq-3a-1}, \varepsilonqref{eq-6-1}, \varepsilonqref{eq-7-1}
and on DP equations \varepsilonqref{dualdpesimp1bb}, \varepsilonqref{dualdpesimp2bb}, \varepsilonqref{eq-7simpbb} for $\gamma=1$, $100$, and $1000$,
on  an instance of the inventory problem with $T=20$ and $N_t=20$.
The corresponding graphs of  $V_2$ (bold dark solid line) and of $V_2^{\gamma}$ for $\gamma=1, 10$, $1000$,
are represented in Figure \ref{fig1inventbisb}. We observe that all functions $V_2^{\gamma}$ are, as expected,
concave upper bounding functions for $V_2$ finite everywhere. We also see that on the domain of
$V_2$, $V_2^{\gamma}$ gets closer to $V_2$ when $\gamma$ increases and eventually coincides with $V_2$
on this domain when $\gamma$ is sufficiently large. Similar graphs were observed for remaining functions
$V_t, V_t^{\gamma}$, $t=3,\ldots,T$.
 \begin{equation}gin{figure}
 \centering
 \begin{equation}gin{tabular}{c}
 \includegraphics[scale=0.5]{BellmanT2.pdf}
 \varepsilonnd{tabular}
\caption{Graph of $V_2$ and of $V_2^{\gamma}$ for
$\gamma=1, 100, 1000$.}\label{fig1inventbisb}
\varepsilonnd{figure}
Therefore, convergence of Dual SDDP with penalizations requires
the coefficients $v_{t k}$ to become arbitrarily large. Proof of the following theorem is given in the Appendix.
\begin{equation}gin{theorem}\label{thconvdualsddppenal}
Consider optimization problem \varepsilonqref{eq-1} and Dual SDDP
with penalizations applied to the dual of this problem.
Let Assumptions (A1) and (A2) hold.
Assume that samples
${ x}i_t^{\varepsilonll}$, $t=2,\ldots,T$, $\varepsilonll \geq 1$,
 in the forward passes are independent,
that $v_{t k+1}\geq v_{t k}$ for all $t,k$,
 and that
$\lim_{k \rightarrow +\infty}
v_{t k} = +\infty$ for all stage $t$. Then the sequence
$V^{k}$ is a deterministic sequence of upper bounds on the optimal value
of \varepsilonqref{eq-1} which converges almost surely to the optimal value of this problem.
\varepsilonnd{theorem}
\par {\textbf{Dual SDDP with feasibility cuts.}}
For dual problems not satisfying the RCR assumption, a subproblem for a given stage $t$
in the forward pass can be infeasible. In this situation, as was done in
Section 5 of \cite{guiguessiopt2016} for SDDP, we can build a feasibility cut
for stage $t-1$ and go back to the previous stage $t-1$ to resolve the problem
with that feasibility cut added, and so on until a sequence of feasible states
is obtained for all stages. In this context, no penalized slack
variables are used, neither in the forward nor in the backward pass.
Since the adaptations from \cite{guiguessiopt2016}
are simple, we skip the details of the derivations
of this SDDP method applied to the dual.
It will be tested in the numerical experiments of Section
\ref{numer}.

\if{

More precisely, we define for stage $t$ the set
$
\mathcal{S}_t:=\{ \pi_t: {\tilde C}_t^\top \pi_t   \leq {\tilde d}_t \}$
where matrices ${\tilde C}_t$ and vectors $\tilde d_t$ are updated along the iterations and such that
a feasible $\pi_t$ must belong to $\mathcal{S}_t$.
The detailed steps of the algorithm are given below.\\

\par {\textbf{Initialization.}} For $t=2,\ldots,T,$ take for $V_t^0$ an affine upper bounding function
for $V_t$. Set iteration counter $k$ to 1.\\
\par {\textbf{Step 1: forward pass (computation of dual trial points)}}.
\par We set $t=1$.
\par {\textbf{While }} $t\leq T$
\par $\hspace*{0.6cm}${\textbf{if }}$t=1$ we compute an optimal solution $\pi_1^k$ of
$$
V^{k-1} =
\begin{equation}gin{array}{l}
\displaystyle \max_{\pi_1} \; b_1^\top \pi_1 + V_{2}^{k-1}(\pi_1)\\
{\underline{\pi}}_1 \leq \pi_1 \leq {\overline{\pi}}_1,
\varepsilonnd{array}
$$
\par $\hspace*{1.2cm}$where $\underline{\pi}_1$, $\overline{\pi}_1$ are as in Lemma \ref{boundedpi} and do $t \leftarrow t+1$.
\par $\hspace*{0.6cm}${\textbf{else }}

\par $\hspace*{1.2cm}${\textbf{if }}$t=T$
\par $\hspace*{1.4cm}$//Defining function $\mathcal{V}_{T}$ by
\begin{equation}gin{equation}\label{forwardfeas0}
\mathcal{V}_{T}(\pi_{T-1}):=
\begin{equation}gin{array}{l}
  \min \limits_{\pi_{T 1},\ldots,\pi_{T N_T}, { z}eta_T} \;
{\textbf{e}}^\top { z}eta_T \\
A_T^{\top} \pi_{T j} \leq c_T,\;j=1,\ldots,N_T,\\
{ z}eta_T \geq 0,
A_{T-1}^\top \pi_{T-1} +
\mathcal{S}um_{j=1}^{N_T} p_{T j} B_{T j}^{\top} \pi_{T j}
\leq c_{T-1}+{ z}eta_T,
  \varepsilonnd{array}
\varepsilonnd{equation}
\par $\hspace*{1.4cm}$//where as before $\textbf{e}$ is a vector of ones then if $\mathcal{V}_{T}(\pi_{T-1}^k)>0$ we have
\par $\hspace*{1.4cm}$//that $\pi_{T-1}^k$ is not feasible
for stage $T-1$.
\par $\hspace*{1.4cm}$//For any $\pi_{T-1}$, problem \varepsilonqref{forwardfeas0} is feasible with finite optimal value
\par $\hspace*{1.4cm}$//$\mathcal{V}_{T}(\pi_{T-1})$
which can be expressed as the optimal value of the
\par $\hspace*{1.4cm}$//corresponding dual problem, i.e.,
\begin{equation}gin{equation}\label{forwardfeas1}
\mathcal{V}_{T}(\pi_{T-1}):=
\begin{equation}gin{array}{l}
  \max \limits_{\delta, \alpha} \; -\mathcal{S}um\limits_{j=1}^{N_{T}}
 c_T^{\top} \alpha_j + \delta^\top (-c_{T-1}+A_{T-1}^{\top} \pi_{T-1})\\
\alpha_j \geq 0, p_{T j} B_{T j} \delta + A_T \alpha_j  = 0,\;j=1,\ldots,N_T,\\
0 \leq \delta \leq e.
  \varepsilonnd{array}
\varepsilonnd{equation}
\par $\hspace*{1.8cm}$Let $\delta,\alpha$ be an optimal
solution of \varepsilonqref{forwardfeas1}
written for $\pi_{T-1}=\pi_{T-1}^k$
\par $\hspace*{1.8cm}$with optimal value $\mathcal{V}_{T}(\pi_{T-1}^k)$.
\par $\hspace*{1.8cm}${\textbf{if }}$\mathcal{V}_{T}(\pi_{T-1}^k)>0$ then we
add a feasibility cut adding
to ${\tilde C}_{T-1}^{\top}$ the
\par $\hspace*{1.8cm}$row $(A_{T-1} \delta)^{\top}$ and to $\tilde d_{T-1}$
the component $\delta^\top c_{T-1} +
c_T^\top \mathcal{S}um\limits_{j=1}^{N_{t}} \alpha_j$.
We
\par $\hspace*{1.8cm}$then go back to the previous stage setting $t \leftarrow t-1$. The corresponding
\par $\hspace*{1.8cm}$feasibility cut
is valid since for $\pi_{T-1}$ feasible for stage
$T-1$ we have
$$\mathcal{V}_{T}(\pi_{T-1})=0 \geq -\mathcal{S}um\limits_{j=1}^{N_{T}}
 c_T^{\top} \alpha_j + \delta^\top (-c_{T-1}+A_{T-1}^{\top} \pi_{T-1}).$$

\par $\hspace*{1.8cm}${\textbf{else }}$t \leftarrow t+1$
\par $\hspace*{1.8cm}${\textbf{end if}}
\par $\hspace*{1.2cm}${\textbf{else }}
\par $\hspace*{1.8cm}$//Let $\mathcal{V}_t$ be the function given by
\begin{equation}gin{equation}\label{forwardfeas2p}
\mathcal{V}_t(\pi_{t-1}):=\begin{equation}gin{array}{l}
  \min \limits_{{ z}eta_t \geq 0,\pi_{t 1},\ldots,\pi_{t N_t}} \; {\textbf{e}}^{\top} { z}eta_t\\
A_{t-1}^{\top} \pi_{t-1}
+ \mathcal{S}um_{j=1}^{N_t} p_{t j} B_{t j}^{\top} \pi_{t j}
\leq c_{t-1}+ { z}eta_t,\\
{\tilde C}_t^{\top} \pi_{t j} \leq  \tilde d_t, j=1,\ldots,N_t.
  \varepsilonnd{array}
\varepsilonnd{equation}
\par $\hspace*{1.8cm}$//Linear program \varepsilonqref{forwardfeas2p} is feasible
and has a finite optimal value for
\par $\hspace*{1.8cm}$//every $\pi_{t-1}$. Therefore
by linear programming duality we can express
\par $\hspace*{1.8cm}$//$\mathcal{V}_t(\pi_{t-1})$
as the optimal value of the corresponding dual problem \par $\hspace*{1.8cm}$//given by
\begin{equation}gin{equation}\label{forwardfeas2}
\mathcal{V}_t(\pi_{t-1}):=\begin{equation}gin{array}{l}
  \max \limits_{\delta, \lambda} \; -\mathcal{S}um\limits_{j=1}^{N_{t}}
 {\tilde d}_t^{\top} \lambda_j +
 \delta^\top (-c_{t-1} + A_{t-1}^{\top} \pi_{t-1})\\
p_{t j} B_{t j} \delta + {\tilde C}_t {\lambda}_j = 0,j=1,\ldots,N_t,\\
0 \leq \delta \leq {\textbf{e}}, \lambda_j \geq 0,\;j=1,\ldots,N_t.
  \varepsilonnd{array}
\varepsilonnd{equation}
\par $\hspace*{1.8cm}$We compute an optimal solution
$\delta, \lambda$ of \varepsilonqref{forwardfeas2} written
with
\par $\hspace*{1.8cm}$$\pi_{t-1}=\pi_{t-1}^k$ with optimal value
$\mathcal{V}_t(\pi_{t-1}^k)$.
\par $\hspace*{1.8cm}${\textbf{if }}$\mathcal{V}_{t}(\pi_{t-1}^k)>0$ then we add
add to ${\tilde C}_{t-1}^{\top}$ the row $(A_{t-1} \delta)^{\top}$, we add
\par $\hspace*{1.8cm}$to $\tilde d_{t-1}$
the component $\delta^\top c_{t-1}+\mathcal{S}um\limits_{j=1}^{N_{t}} {\tilde d}_t^{\top} \lambda_j$, and do $t \leftarrow t-1$. The
\par $\hspace*{1.8cm}$corresponding
feasibility cut
is valid since for $\pi_{t-1}$ feasible for stage
\par $\hspace*{1.8cm}$$t-1$ we have
$$\mathcal{V}_{t}(\pi_{t-1})=0 \geq -\mathcal{S}um\limits_{j=1}^{N_{t}}
 {\tilde d}_t^{\top} \lambda_j +
 \delta^\top (-c_{t-1} + A_{t-1}^{\top} \pi_{t-1}).$$
\par $\hspace*{1.8cm}${\textbf{else }}we compute an optimal solution $\pi_{t}^k$ of
\begin{equation}gin{equation}\label{forwardfeas3}
\begin{equation}gin{array}{cll}
  \max\limits_{\pi_{t 1},\ldots,\pi_{t N_t}} &  \mathcal{S}um\limits_{j=1}^{N_{t}}  p_{t j}  \left[ b_{t j}^\top
  \pi_{t j}+ {V}_{t+1}^{k-1}(\pi_{t j})\right]\\
    {\rm s.t.}&
  A_{t-1}^\top \pi_{t-1}^k +
  \mathcal{S}um\limits_{j=1}^{N_{t}} p_{tj} B_{t j}^\top  \pi_{t j}  \le c_{t-1},\\
  & {\underline{\pi}}_t \leq \pi_{t j} \leq {\overline{\pi}}_t, {\tilde C}_t^{\top} \pi_{t j} \leq \tilde d_t,\;j=1,\ldots,N_t.
   \varepsilonnd{array}
\varepsilonnd{equation}
\par $\hspace*{1.8cm}$An optimal solution of the problem above has $N_t$ components
$(\pi_{t 1}, \pi_{t 2}$,
\par $\hspace*{1.8cm}$ $\ldots,\pi_{t N_t})$
for $\pi_t$. We generate a realization $\tilde { x}i_t^k$ of
${ x}i_t^k \mathcal{S}im { x}i_t$
\par $\hspace*{1.8cm}$independently of previous realizations
$\tilde { x}i_2^1,\ldots$, $\tilde { x}i_{T-1}^1$,$\ldots$,
$\tilde { x}i_2^k,\ldots,\tilde { x}i_{t-1}^k$,
\par $\hspace*{1.8cm}$and take $\pi_{t}^k = \pi_{t j_t(k)}$ where index $j_t(k)$ satisfies
$\tilde { x}i_{t}^k = { x}i_{t j_t(k)}$.
\par $\hspace*{1.8cm}$We then go to the next
stage: $t \leftarrow t+1$.
\par $\hspace*{1.8cm}${\textbf{end if}}
\par $\hspace*{1.2cm}${\textbf{end if}}
\par $\hspace*{0.6cm}${\textbf{end if}}
\par {\textbf{end while}}\\
\par {\textbf{Step 2: backward pass}}. For
$t=2,\ldots,T$, we compute a new cut for $V_t$ to build $V_t^k$ exactly as in Dual SDDP
with penalizations taking penalizations
$v_{t k}=0$ (which is possible since trial points
are feasible).\\
\par {\textbf{Step 3:}} Do $k \leftarrow k+1$ and go to Step 1.\\

In Dual SDDP with feasibility cuts
all subproblems solved in the backward passes are feasible and $V^k$ is a nonincreasing deterministic
 sequence of upper bounds on the optimal value of \varepsilonqref{eq-1}.
Moreover, following the proof of
Theorem 5.1 in \cite{guiguessiopt2016}, we obtain
that if the samples in the forward passes
are independent then
the sequence
$V^{k}$ converges to the optimal value of
primal problem \varepsilonqref{eq-1}.

}\fi

\mathcal{S}ubsection{Dual SDDP for problems with uncertainty in all parameters}

We have seen in Section \ref{wrdpequ} how to write DP equations
on the dual problem of a MSLP when all data $(A_t,B_t,c_t,b_t)$ in $({ x}i_t)$
is random. In this situation, cost-to-go functions $V_t$ are functions
$V_t(\pi_{t-1},{ x}i_{t-1})$ of both past decision $\pi_{t-1}$
and past value ${ x}i_{t-1}$ of process $({ x}i_t)$.
Also recall that functions
$V_t(\cdot,{ x}i_{t-1})$ are concave for all ${ x}i_{t-1}$.
Therefore, Dual SDDP with penalizations from the previous section
must be modified as follows. For each stage $t=2,\ldots,T,$ instead of computing
just one approximation of a single function (function $V_t$), we now need
to compute approximations of $N_t$ functions, namely
concave cost-to-go functions $V_t(\cdot,{ x}i_{t-1 j})$, $j=1,\ldots,N_t$.
The approximation $V_{t j}^k$ computed for $V_t(\cdot,{ x}i_{t-1 j})$  at iteration $k$ is a polyhedral function
$V_{t j}^k$ given by:
$$
V_{t j}^k( \pi_{t-1}) = \displaystyle \min_{0 \leq i \leq k}  \;{\overline{\theta}}_{t j}^i + \langlegle {\overline{\begin{equation}ta}}_{t j}^i , \pi_{t-1} {\ranglegle}ngle.
$$
Therefore more computational effort is needed.
However, the adaptations of the method can be easily written.
More specifically, at iteration $k$, in the forward pass, dual trial points
are obtained replacing $V_{t}(\cdot,{ x}i_{t-1 j})$ by $V_{t j}^{k-1}$
and in the backward pass a cut is computed at stage $t$ for $V_{t}(\cdot,{ x}i_{t-1 j_k})$ with
$j_k$ satisfying ${ x}i_{t-1 j_k}=\tilde { x}i_{t-1}^k$ where
$\tilde { x}i_{t-1}^k$ is the sampled value of ${ x}i_{t-1}$ at iteration $k$.

\mathcal{S}ubsection{Dual SDDP for problems with interstage dependent cost coefficients}
\label{dualsddp}

We consider problems of form
\varepsilonqref{eq-1} where costs  $c_t$ affinely depend on their  past
while  $b_t$  are stagewise  independent.
Specifically, similar to derivations of Section \ref{sec-invent},
suppose that $c_t$ follow  a multiplicative vector autoregressive process of form
\begin{equation}gin{equation}\label{inv-2}
\begin{equation}gin{array}{ll}
 c_t = {\rm Var}epsilon_t \circ \left(  \mathcal{S}um_{j=1}^p \Phi_{t j} c_{t-j} + \mu_t  \right),
 \varepsilonnd{array}
\varepsilonnd{equation}
with  $(x \circ y)_i=x_i y_i$ denoting the componentwise product,  and where matrices $\Phi_{t j}$ and vectors  $\mu_t \geq 0$ as well as $c_1,\ldots,c_{2-p} \geq 0$ are given.

We assume that the  process
$(b_t,{\rm Var}epsilon_t)$ is stagewise
independent and that the support
of $b_t, {\rm Var}epsilon_t$ is the finite
set
$$\{
(b_{t 1},{\rm Var}epsilon_{t 1}),
\ldots,
(b_{t N_t},{\rm Var}epsilon_{t N_t})
\},$$
with
${\rm Var}epsilon_{t i}>0$
and  $p_{t i}=
\mathbb{P}\{
(b_t, {\rm Var}epsilon_t)=
(b_{t i}, {\rm Var}epsilon_{t i})\}$,
$i=1,\ldots,N_t$.
For some values of $\Phi_{t j}$ (for instance for matrices with
nonnegative entries), this guarantees that all realizations of the price process $\{c_t\}$ are positive.
The developments which follow
can be easily extended to other linear models for $\{c_t\}$, for instance SARIMA or PAR
models, see \cite{vgcoap2014} for the definition of state vectors of minimal
size for generalized linear models.

Using the notation
$c_{t_1:t_2}=(c_{t_1},c_{t_1+1},\ldots,c_{t_2-1},c_{t_2})$ for
$t_1 \leq t_2$ integer, for the corresponding primal problem (of the form \varepsilonqref{eq-1}),  we can write the following Dynamic Programming equations:
define $\mathcal{Q}_{T+1} \varepsilonquiv 0$ and for $t=2,\ldots,T$,
\begin{equation}gin{equation} \label{defqtb}
\mathcal{Q}_t ( x_{t-1} , c_{t-p:t-1} ) =
\mathbb{E}_{b_t, {\rm Var}epsilon_t}{\cal B}ig[
Q_t ( x_{t-1} , c_{t-p:t-1}, b_t, {\rm Var}epsilon_t ) {\cal B}ig]
\varepsilonnd{equation}
where $Q_t (x_{t-1} , c_{t-p:t-1}, b_t, {\rm Var}epsilon_t )$ is given by
\begin{equation}gin{equation}\label{deffraktb}
\begin{equation}gin{array}{ll}
\displaystyle \min_{x_t \geq 0}  {\cal B}ig[{\rm Var}epsilon_t \circ {\cal B}ig(\displaystyle \mathcal{S}um_{j=1}^p \Phi_{t j} c_{t-j} + \mu_t {\cal B}ig) {\cal B}ig]^{\top}
x_t + \mathcal{Q}_{t+1}{\cal B}ig(x_t,
c_{t+1-p:t-1},{\rm Var}epsilon_t \circ {\cal B}ig(\displaystyle \mathcal{S}um_{j=1}^p \Phi_{t j} c_{t-j} + \mu_t {\cal B}ig){\cal B}ig ) \\
 A_t x_t + B_t x_{t-1} = b_t,
\varepsilonnd{array}
\varepsilonnd{equation}
while the first stage problem is
$$
\begin{equation}gin{array}{l}
\displaystyle \min_{x_1 \geq 0} \; c_1^\top x_1 + \mathcal{Q}_2( x_1,c_{2-p:1}) \\
A_1 x_1 = b_1.
\varepsilonnd{array}
$$

Standard  SDDP does not apply directly to solve Dynamic Programming equations
\varepsilonqref{defqtb}-\varepsilonqref{deffraktb} because functions $\mathcal{Q}_t$ given by \varepsilonqref{defqtb}-\varepsilonqref{deffraktb} are not convex.
Nevertheless,   we can use the Markov Chain discretization variant of SDDP to solve
Dynamic Programming equations \varepsilonqref{defqtb}-\varepsilonqref{deffraktb}.
On the other hand, as pointed above,  it is possible to apply SDDP for the dual problem with the added state variables.
Along the lines of Section \ref{wrdpequ}
we can write Dynamic Programming equations for the dual,
now with function $V_t$ depending on
$\pi_{t-1},c_{t-1},\ldots,c_{t-p}$.

\if{
The following Dynamic Programming equations
for the dual of \varepsilonqref{eq-1} with $(c_t)$ of form \varepsilonqref{inv-2} can be be written.
For the last stage $T$, we have to solve the problem:
\begin{equation}gin{equation}\label{dualdpesimp1b}
\begin{equation}gin{array}{l}
  \max\limits_{\pi_{T 1},\ldots,\pi_{T N_T}} \; \mathcal{S}um\limits_{j=1}^{N_T} p_{T j} \pi_{T j}^{\top} b_{T j}\\
    A_{T}^{\top} \pi_{T j} \leq
    {\rm Var}epsilon_{T j} \circ \left(\mu_{T} +  \mathcal{S}um_{\varepsilonll=1}^p  \Phi_{T \varepsilonll} c_{T-\varepsilonll} \right),\;j=1,...,N_T,\\
      \mathcal{S}um_{j=1}^{N_T} p_{T j} B_{T}^{\top} \pi_{T j} \leq c_{T-1}-A_{T-1}^{\top} \pi_{T-1},
   \varepsilonnd{array}
\varepsilonnd{equation}
with optimal value $V_T(\pi_{T-1} , c_{T-p:T-1})$.

Next for stage $t=2,\ldots,T-1$, given $V_{t+1}$, we need
to solve the problem
\begin{equation}gin{equation}\label{dualdpesimp2b}
\begin{equation}gin{array}{cl}
 \max \limits_{\pi_{t 1},\ldots,\pi_{t N_t}} & \mathcal{S}um \limits_{j=1}^{N_{t}}  p_{t j}
\left( \pi_{t j}^\top b_{t j} + V_{t+1}{\cal B}ig(\pi_{t j},  c_{t+1-p:t-1},{\rm Var}epsilon_{t j} \circ (
\mu_{t} +  \mathcal{S}um_{\varepsilonll=1}^p  \Phi_{t \varepsilonll} c_{t-\varepsilonll} ){\cal B}ig) \right)\\
{\rm s.t.} &  \mathcal{S}um_{j=1}^{N_{t}} p_{t j}
  B_{t}^{\top} \pi_{t j}  \leq c_{t-1} - A_{t-1}^{\top} \pi_{t-1},
   \varepsilonnd{array}
   \varepsilonnd{equation}
while the first stage problem is
 \begin{equation}gin{equation}\label{eq-7simpb}
  \max_{\pi_1} \;\;      \pi_1^\top  b_1  +  V_{2}(\pi_1 , c_{2-p:1}).
\varepsilonnd{equation}
}\fi
These functions are concave
and  therefore we can apply
Dual SDDP with penalizations to these DP equations
to build polyhedral approximations of these functions $V_t$ of form
\begin{equation}gin{equation}
V_t^k(\pi_{t-1},c_{t-1},\ldots,c_{t-p})
= \min_{0 \leq i \leq k} \theta_{t}^i +  \langlegle \begin{equation}ta_{t 0}^i , \pi_{t-1} {\ranglegle}ngle + \mathcal{S}um_{j=1}^p  \langlegle \begin{equation}ta_{t j}^i , c_{t-j} {\ranglegle}ngle
\varepsilonnd{equation}
at iteration $k$.

\mathcal{S}etcounter{equation}{0}
\mathcal{S}ection{Numerical experiments}
\label{numer}

In this section, we report numerical results obtained applying Primal SDDP and  variants of Dual SDDP

to the  inventory problem and to the  Brazilian interconnected power system problem.
All methods were implemented in Matlab and run on an Intel Core i7, 1.8GHz, processor
with 12,0 Go of RAM. Optimization problems were solved using Mosek \cite{mosek}.
 \if{
When the data process  of problem \varepsilonqref{eq-1} is  interstage dependent, in some cases the problem can be reformulated  in a stagewise independent form by adding  appropriate state variables. In case the data process can   represented as an autoregressive process and only the  right hand sides $b_t$  of \varepsilonqref{eq-1} are random, this preserves the linearity of the obtained problem (see  section \ref{sec-invent}). On the other hand if the costs $c_t$ are given by an
autoregressive process, then such approach of adding state variables destroys the linearity and convexity of the primal problem. It could be noted that  in the dual problem \varepsilonqref{eq-2} the costs are in the right hand side  of the constraints,  and   hence such approach can be applied to the dual problem.
}\fi

\mathcal{S}ubsection{Dual SDDP for the inventory problem}\label{sec:inventorysim}

\if{
namely as
$$
-a_1 x_0 +
\left\{
\begin{equation}gin{array}{l}
\displaystyle \min \mathbb{E}_{{ x}i}\left[  \mathcal{S}um_{t=1}^T c_t^{\top} x_t + d_t^{\top} u_t  \right] \\
A_1 x_1 + C_1 u_1 = b_1, D_1 x_1 + F_1 u_1 \leq  f_1,\\
A_t x_t + B_t x_{t-1} + C_t u_t = b_t, \mbox{a.s.},\;t=2,\ldots,T,\\
D_t x_t + E_t x_{t-1} + F_t u_t \leq  f_t,\mbox{a.s.},\;t=2,\ldots,T,
\varepsilonnd{array}
\right.
$$
where the underlying stochastic process $\mathcal{D}_t$ is the demand process, control $u_t=[y_t;w_t;v_t]$,
$$
\begin{equation}gin{array}{lllll}
c_t=  \left\{
\begin{equation}gin{array}{ll}
-a_{t+1} & \mbox{if }t<T,\\
0 & \mbox{for }t=T,
\varepsilonnd{array}
\right.
&
A_1= [-1],
&
C_1=[1 \;  0  \; 0],
&
D_1=[0;0;0;0;0].
&
F_1=\left[
\begin{equation}gin{array}{ccc}
0 & -1 & 0 \\
-1 & -1 & 0 \\
0 & 0 & -1 \\
1 & 0 & -1 \\
-1 & 0 & 0
\varepsilonnd{array}
\right],
\varepsilonnd{array}
$$
$$
\begin{equation}gin{array}{llll}
d_t=\left[
\begin{equation}gin{array}{c}
a_t \\
g_t \\
h_t
\varepsilonnd{array}
\right],
b_t=[\mathcal{D}_t], &
f_1= \left[
\begin{equation}gin{array}{c}
0\\
-{ x}i_1\\
0\\
{ x}i_1 \\
-x_0
\varepsilonnd{array}
\right],
 &   f_t= \left[
\begin{equation}gin{array}{c}
0\\
-\mathcal{D}_t\\
0\\
\mathcal{D}_t \\
0
\varepsilonnd{array}
\right], t \geq 2,
\varepsilonnd{array}
$$
and for $t=2,\ldots,T$:
$$
\begin{equation}gin{array}{llllll}
A_t=[-1],&B_t=[0],&C_t=[1 \; 0 \; 0],&D_t=[0;0;0;0;0]&E_t=[0;0;0;0;1],&F_t=F_1.
\varepsilonnd{array}
$$
}\fi
We consider the inventory problem \varepsilonqref{inv-1}
with parameters $a_t=1.5+\cos(\frac{\pi t}{6})$, $p_{t i}=\frac{1}{N}$ where $N$ is the number of realizations for each stage, ${ x}i_{t j}=(5+0.5 t)(1.5+0.1 z_{t j})$ where
$(z_{t 1},\ldots,z_{t N})$ is a sample from the standard Gaussian
distribution, $x_0=10$, $g_t=2.8$, and
$h_t=0.2$.

\par {\textbf{Illustrating the correctness of DP equations \varepsilonqref{eq-3a-1}, \varepsilonqref{eq-6-1}, \varepsilonqref{eq-7-1} and checking the convergence of the variants of Dual SDDP.}}
\begin{equation}gin{figure}
 \centering
 \begin{equation}gin{tabular}{cc}

  \includegraphics[scale=0.48]{CutsBellmanT2100} &
\includegraphics[scale=0.48]{CutsBellmanT2FCuts}
 \varepsilonnd{tabular}
 \caption{Graph of $V_2$ (bold black solid line) and cuts computed for $V_2$
by Dual SDDP with penalizations $v_{t k}=100$ (left panel) and
Dual SDDP with feasibility cuts (right panel).}\label{fig1invent}
\varepsilonnd{figure}
We solve this inventory problem using Dynamic Programming applied both
to DP equations \varepsilonqref{eq-3a-1}, \varepsilonqref{eq-6-1}, \varepsilonqref{eq-7-1} and to DP equations \varepsilonqref{dualdpesimp1bb}-\varepsilonqref{dualdpesimp2bb}
for $\gamma=1, 10, 1000$.
In this latter case, we obtain approximations of functions $V_t^{\gamma}$.
We also run Primal SDDP, Dual SDDP with feasibility cuts,
and Dual SDDP with penalties $v_{t k}=1$, $10$, $1000$,
on the same instance, knowing that Dual SDDP variants were run
for 100 iterations (the upper bounds computed by these methods
stabilize in less than 10 iterations) and Primal SDDP was stopped
when
the gap is $<0.1$ where the gap is defined as $\frac{Ub-Lb}{Ub}$ where $Ub$ and $Lb$ correspond to upper and lower bounds
 computed
by Primal SDDP along iterations.
The lower bound $Lb$ is the optimal value of the first stage problem and the
upper  bound $Ub$ is the upper end of a 97.5\%-one-sided confidence interval on the optimal value obtained using the sample
of total costs computed by all previous forward passes. With this stopping criterion and the considered instance
of the inventory problem, Primal SDDP was run for 232 iterations.

In Figure \ref{fig1invent}, we report the graph of $V_2$ and the cuts computed for
$V_2$ by Dual SDDP with feasibility cuts (right panel) and
Dual SDDP with penalties $v_{t k}=100$ (left panel).
All cuts are, as expected, upper bounding affine functions for $V_2$
on its domain. However, it is interesting to notice that for Dual SDDP with feasibility
cuts, few different cuts are computed and these cuts are tangent or very close to $V_2$
at the trial points.
On the contrary, Dual SDDP with penalties
may compute many cuts dominated by others on the domain of $V_2$.
Therefore, cut selection techniques, for instance along the lines of
\cite{guiguesejor17} \cite{guiguesbandarra18} using Limited Memory Level 1 cut selection, could be useful for Dual SDDP.

We report in Table \ref{table0}
the approximate optimal values  and the time needed to compute them
with Primal SDDP, Dual SDDP, and Dynamic Programming
applied to respectively \varepsilonqref{eq-3a-1}, \varepsilonqref{eq-6-1}, \varepsilonqref{eq-7-1}
and \varepsilonqref{dualdpesimp1bb}, \varepsilonqref{dualdpesimp2bb}, \varepsilonqref{eq-7simpbb} with $\gamma=1, 100, 1000$.
The approximate optimal values reported are the last upper bound computed for variants of Dual SDDP
and the last lower bound computed for Primal SDDP. All approximate optimal values
are very close (showing that all variants were correctly implemented) and Dynamic Programming is much slower than the other sampling-based algorithms.
For Dual SDDP with penalization, if penalties are too small the upper bound can be $+\infty$ while if penalties
are sufficiently large the algorithm converges to an optimal policy.

\begin{equation}gin{table}
\centering
\begin{equation}gin{tabular}{|c|c|c|}
\hline
{\tt{Method}} & {\tt{Optimal value}} &  {\tt{CPU time (s.)}}    \\
\hline
{\tt{DP}} on \varepsilonqref{eq-3a-1}, \varepsilonqref{eq-6-1}, \varepsilonqref{eq-7-1}  &   321.6& 685\\
\hline
{\tt{DP}} on  \varepsilonqref{dualdpesimp1bb}, \varepsilonqref{dualdpesimp2bb}, \varepsilonqref{eq-7simpbb}, $\gamma=1$  & $+\infty$ & 2 860 \\
\hline
{\tt{DP}} on  \varepsilonqref{dualdpesimp1bb}, \varepsilonqref{dualdpesimp2bb}, \varepsilonqref{eq-7simpbb}, $\gamma=100$  &  322.2 & 3 808\\
\hline
{\tt{DP}} on \varepsilonqref{dualdpesimp1bb}, \varepsilonqref{dualdpesimp2bb}, \varepsilonqref{eq-7simpbb}, $\gamma=1000$  & 321.8 & 3 376 \\
\hline
{\tt{Primal SDDP}}   &  322.5  &     105     \\
\hline
{\tt{Dual SDDP with penalties}}, $v_{t k}=1$  & 2 131.4 & 9.4\\
\hline
{\tt{Dual SDDP with penalties}}, $v_{t k}=100$  & 322.5& 11.3\\
\hline
{\tt{Dual SDDP with penalties}}, $v_{t k}=1000$  &  322.5 & 11.9\\
\hline
{\tt{Dual SDDP with feasibility cuts}}  &  322.5 &10.6 \\
\hline
\varepsilonnd{tabular}
\caption{Optimal value and CPU time needed (in seconds) to compute them on an instance of the inventory problem
with $T=N_t=20$ by Dynamic Programming (DP), Primal SDDP, and variants of Dual SDDP.}\label{table0}
\varepsilonnd{table}
Finally, we report for this instance in Figure \ref{fig1bounds}
the evolution of the lower bound $Lb$ and upper bound $Ub$ computed by Primal SDDP and the upper bounds computed by
Dual SDDP with penalties $v_{t k}=1000$ and Dual SDDP with feasibility cuts.
With  Dual SDDP, the upper bound is naturally large at the first iteration but decreases much quicker than the upper
bound $Ub$ computed by Primal SDDP, especially for Dual SDDP with feasibility cuts, with all upper bounds
converging to the optimal value of the problem.\\

\begin{equation}gin{figure}
 \centering
 \begin{equation}gin{tabular}{cc}
\includegraphics[scale=0.48]{BoundsT20N20First.pdf}
 &
\includegraphics[scale=0.48]{BoundsT20N20Second.pdf}\\
 \varepsilonnd{tabular}
\caption{Left: upper and lower bounds computed by Primal SDDP
and upper bounds computed by Dual SDDP with feasibility cuts
and Dual SDDP with penalties $v_{t k}=1000$
for the first 10 iterations.
Right: same outputs for iterations 10,$\ldots,100$.}\label{fig1bounds}
\varepsilonnd{figure}

\par {\textbf{Tests on a larger instance.}} We now run Primal and Dual SDDP on a larger instance
with $T=100$ and $N_t=100$ for 600 iterations.
The evolution of the upper bounds computed along the iterations
of Dual SDDP (both with feasibility cuts and with penalizations $v_{t k}=1000$)
and of the upper and lower bounds computed by Primal SDDP
are reported in Table \ref{tablebounditer}
for
iterations $2, 3, 5, 10, 50, 100, 200, 300, 400, 500,$ and $600$.
We see that for the first iterations, the upper bound decreases more quickly with the variants
of Dual SDDP, the most important decrease being obtained for Dual SDDP with feasibility cuts.
However, on this instance, the convergence of Dual SDDP with feasibility cuts is slower, i.e.,
a solution of high accuracy is obtained quicker using Dual SDDP with  penalizations.
More precisely, we fix confidence levels ${\rm Var}epsilon=0.2, 0.15, 0.1, 0.05, 0.01$, and
for each confidence level, we compute the time needed, running Primal and Dual SDDP in parallel, to obtain a solution with relative accuracy ${\rm Var}epsilon$ stopping the algorithm when
the upper bound {\tt{Ub\_D}} computed by a variant of Dual SDDP and the lower bound Lb, computed by Primal SDDP,
satisfies ({\tt{Ub\_D}}-Lb)/{\tt{Ub\_D}}$<{\rm Var}epsilon$. The results are reported in Table \ref{table1}.
In this table, we also report the time needed to obtain a solution of relative accuracy
${\rm Var}epsilon$ using only the information provided by Primal SDDP, stopping the algorithm when
(Ub-Lb)/Ub$<{\rm Var}epsilon$.

We observe that if ${\rm Var}epsilon$ is not too small, the smallest CPU time
is obtained combining Primal SDDP with Dual SDDP
with feasibility cuts while when ${\rm Var}epsilon$ is small (0.05 and 0.01) the smallest CPU
time is obtained combining Primal SDDP with Dual SDDP with penalizations. For ${\rm Var}epsilon=0.05$ and $0.01$,
600 iterations are even not enough to get a solution of relative accuracy ${\rm Var}epsilon$ using
Primal SDDP or combining Primal SDDP and Dual SDDP with feasibility cuts.

\if{
\begin{equation}gin{figure}
{\color{blue} this figure takes the whole page, can we remove it?}
 \centering
 \begin{equation}gin{tabular}{cc}
\includegraphics[scale=0.48]{BoundsT100N100First.pdf}
 &
\includegraphics[scale=0.48]{BoundsT100N100Second.pdf}\\
 \varepsilonnd{tabular}
\caption{Left: upper and lower bounds computed by Primal SDDP
and upper bounds computed by Dual SDDP with feasibility cuts
and Dual SDDP with penalties $v_{t k}=1000$
for the first 10 iterations for an instance of the inventory problem with $T=N_t=100$.
Right: same outputs for iterations 10,$\ldots,100$.}\label{fig2bounds}
\varepsilonnd{figure}
}\fi

\begin{equation}gin{table}[H]
\centering
\begin{equation}gin{tabular}{|c|c|c|c|c|}
\hline
Iteration & \begin{equation}gin{tabular}{c}{\tt{Primal}}\\{\tt{SDDP}}\\Lb\varepsilonnd{tabular} & \begin{equation}gin{tabular}{c}{\tt{Primal}}\\
{\tt{SDDP}}\{\cal U}b\varepsilonnd{tabular}  &  \begin{equation}gin{tabular}{c}{\tt{Dual SDDP with}}\\{\tt{feasibility}}\\{\tt{cuts}}\varepsilonnd{tabular} & \begin{equation}gin{tabular}{c}{\tt{Dual SDDP}}\\{\tt{with}}\\  {\tt{penalties}}\varepsilonnd{tabular}\\
\hline
2  & 656.4 &25 443  & 20 002 &20 015 \\
\hline
3  & 713.1  &19 340  &8 693.1  & 20 012\\
\hline
5  & 3361.8  &  14 800  &7 246.8&19 993\\
\hline
10 & 5330.1   &  10 662 &  5 736.6 & 16 452 \\
\hline
50     &    5483.1    &6 594.5& 5721.8 & 5500.9\\
\hline
100     &  5483.5     & 6 039.2&  5715.1 & 5484.8\\
\hline
200     &  5483.6     &5 762.4 & 5710.0  & 5484.2\\
\hline
300     &   5483.7   & 5 671.0 &  5704.6 & 5484.0\\
\hline
400     &  5483.7    &5 625.3  & 5702.7 & 5483.9\\
\hline
500     &  5483.7    &5 597.9  &  5702.5 & 5483.8\\
\hline
600     &   5483.7   & 5 579.9 &  5702.2 &5483.8 \\
\hline
\varepsilonnd{tabular}
\caption{For an instance of the inventory problem with $T=N_t=100$,
lower bound Lb and upper bound Ub computed by Primal SDDP and upper bounds computed
by Dual SDDP with
feasibility cuts and Dual SDDP with
penalties $v_{t k}=1000$ along iterations.}\label{tablebounditer}
\varepsilonnd{table}

In Figure \ref{figtimeinvent}, we report the cumulative CPU  time along iterations of
 all methods. We see that each iteration requires a similar computational bulk and the CPU time
 increases exponentially with the number of iterations.
 \begin{equation}gin{table}
\centering
\begin{equation}gin{tabular}{|c|c|c|c|}
\hline
${\rm Var}epsilon$ & {\tt{Primal SDDP}} & \begin{equation}gin{tabular}{c}{\tt{Dual SDDP with}}\\{\tt{feasibility cuts}}\varepsilonnd{tabular} & \begin{equation}gin{tabular}{c}{\tt{Dual SDDP with}}\\{\tt{penalties $v_{t k}=1000$}}\varepsilonnd{tabular}\\
\hline
0.2  & 300.2 & 29.5 & 35.8 \\
\hline
0.15  & 459.8  & 35.8 & 41.2 \\
\hline
0.1  & 825.6  & 48.3   &48.3\\
\hline
0.05 & 2366.2   &96.1&  61.5  \\
\hline
0.01     &    -    & -  &103.2 \\
\hline
\varepsilonnd{tabular}\\
\caption{Time needed (in seconds) to obtain a solution of relative accuracy ${\rm Var}epsilon$
with Primal SDDP, Dual SDDP with feasibility cuts, and Dual SDDP with penalties $v_{t k}=1000$
for an instance of the inventory problem with $T=N_t=100$.}\label{table1}
\varepsilonnd{table}

\begin{equation}gin{figure}
 \centering
\includegraphics[scale=0.6]{TimeIter.pdf}
\caption{Cumulative CPU time along iterations of Primal SDDP, Dual SDDP with feasibility cuts, and
Dual SDDP with penalizations $v_{t k}=1000$.}\label{figtimeinvent}
\varepsilonnd{figure}

\if{
Finally, for Dual SDDP with penalizations, we report in Figure \ref{fig3error} the maximal and mean values of
${ z}eta_t^{k}$ (recall that for the inventory problem ${ z}eta_t$ in \varepsilonqref{forwarddualsddppen}
are real-valued)  along iterations for each stage $t$, where ${ z}eta_t^{k}$ is an optimal
value
of ${ z}eta_t$ in \varepsilonqref{forwarddualsddppen} for iteration $k$. The corresponding values
are positive, meaning that indeed RCR does not hold for the dual of the inventory problem.
}\fi

\if{
\begin{equation}gin{figure}
 \centering
 \begin{equation}gin{tabular}{cc}
\includegraphics[scale=0.4]{MaxViolationT100N100.pdf}
 &
\includegraphics[scale=0.4]{MeanViolationT100N100.pdf}\\
 \varepsilonnd{tabular}
\caption{Left: maximal constraint violation as a function of the stage. Right: Mean constraint violation as a function of
the stage. }\label{fig3error}
\varepsilonnd{figure}
}\fi

\mathcal{S}ubsection{Sensitivity analysis for the inventory problem}
\label{invnumer}

Consider the inventory problem of Section \ref{sec:inventorysim} with $(\mathcal{D}_t)$ as in
\varepsilonqref{inv-2-1bis} and $T=10$ stages.
For this problem, the derivatives
from Proposition \ref{pr-valsens}
are given by
\begin{equation}gin{eqnarray}
\label{inv-12}
&\partial {\rm Var}theta(\phi,\mu)/\partial \phi=
\partial L(\mbox{\boldmath$a$}r{x},\mbox{\boldmath$a$}r{y},\mbox{\boldmath$a$}r{\pi})/\partial \phi=
\mbox{\tiny\boldmath$b$}e\left[\mathcal{S}um_{t=1}^T
\mbox{\boldmath$a$}r{ \pi}_t  {\varepsilonpsilon}silon_{t}  \mathcal{D}_{t-1}
 \right] ,\\
 \label{inv-13}
& \partial {\rm Var}theta(\phi,\mu)/\partial \mu=
 \partial L(\mbox{\boldmath$a$}r{x},\mbox{\boldmath$a$}r{y},\mbox{\boldmath$a$}r{\pi})/\partial \mu=
\mbox{\tiny\boldmath$b$}e\left[\mathcal{S}um_{t=1}^T
\mbox{\boldmath$a$}r{ \pi}_t  {\varepsilonpsilon}silon_{t}
 \right] ,
\varepsilonnd{eqnarray}
where $(\mbox{\boldmath$a$}r{x},\mbox{\boldmath$a$}r{y})$ is an optimal solution of the primal  problem and $\mbox{\boldmath$a$}r{\pi}$ are the corresponding Lagrange multipliers.
Our goal is to compute these derivatives
 solving
the primal and dual problems by
respectively Primal and Dual SDDP.

We consider $4$ instances with $(\phi,\mu)=(0.01,0.1)$, $(0.01,3.0)$, $(0.001,0.1)$, and $(0.001,3.0)$.
The remaining parameters of these instances are those from the previous section.
We discretize both the primal and dual problem into $N_t = 100$ samples
for each stage $t = 2,\ldots,10$. We take the relative error
${\rm Var}epsilon = 0.01$ for the stopping criterion and use
$10\,000$ Monte Carlo simulations to estimate the expectations in \varepsilonqref{inv-12}, \varepsilonqref{inv-13}.
For Primal SDDP, the upper bound Ub and lower bound Lb at termination are given
in Table $\ref{tab-sen-num-primal}$ for the four instances.
\begin{equation}gin{table}[htbp]
	\centering
	\mathcal{S}mall
	\begin{equation}gin{tabular}{|c|c|c|c|c|}
		\hline
	Bound		 & {\texttt{Instance 1}} &
			{\texttt{Instance 2}} & {\texttt{Instance 3}} & {\texttt{Instance 4}} \\
			\hline

	   Ub &17.9176    &  478.687 & 15.3940 & 404.242 \\
	   \hline
	    Lb	&17.9163    &475.017  & 15.3927 & 402.913 \\
	    \hline

	\varepsilonnd{tabular}
	\caption{Upper and lower bounds at the last iteration of Primal SDDP.}
	\label{tab-sen-num-primal}
\varepsilonnd{table}

The optimal mean values of Lagrangian multipliers for the demand constraints computed, for a given stage $t \geq 2$, averaging
over the $10\,000$ values obtained simulating $10\,000$ forward passes after termination, are given in Table \ref{tab-sen-num-dual}. In this table, \texttt{LM} stands for the multipliers obtained using Primal SDDP as explained in Remark \ref{remarkcomputedualmpsddp}
whereas \texttt{Dual} stands for the multipliers obtained using Dual SDDP with penalties.
The fact that the multipliers obtained are close for both methods illustrates the validity of the two alternatives we
discussed in Sections \ref{sec-sens}-\ref{sec:dualsddp} to compute derivatives of the value function of a MSP.
\begin{equation}gin{table}[htbp]
	\centering
	\mathcal{S}mall
	\begin{equation}gin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
		Stage & \multicolumn{2}{|c|}{\texttt{Instance 1}} & \multicolumn{2}{|c|}{\texttt{Instance 2}} & \multicolumn{2}{|c|}{\texttt{Instance 3}} & \multicolumn{2}{|c|}{\texttt{Instance 4}} \\
		\hline
		{} &\texttt{LM} & \texttt{Dual} &\texttt{LM} & \texttt{Dual} &\texttt{LM} & \texttt{Dual} &\texttt{LM} & \texttt{Dual}\\
    \hline
        2	 &0.2465
     & 0.2373  &1.6701	&1.66959 & 0.0444	&0.0328 &1.666	&1.666\\
     \hline
        3    &0.3218      &  0.31095 &1.4098	&1.4120 &0.1421	&0.1340 &1.406	&1.409\\
         \hline
        4     & 0.3268   &	0.3221
  &0.9862	&0.9861 &0.19439	&0.18974 &0.984 &0.984\\
   \hline
        5     &0.3086   &	0.3058 & 0.6330	&0.6329 &0.2145	&0.2128  &0.6327	&0.6327\\
         \hline
        6     &0.3408   &	0.3412 &0.49998	& 0.499897   &0.2708	 &0.2717  &0.4999	&0.4998\\
         \hline
        7     &0.5026 	&0.5051   &0.63397	&0.63397   &0.4378	&0.4418 &0.6339	&0.6339\\
         \hline
       	8    &0.7047    & 0.7049
 &0.8348	& 0.8340 &0.6404 &	0.6413 &0.8349	&0.8334\\
  \hline
        9     &0.8985  &	0.9032 &1.0322	&1.0343    &0.83501	&0.8401 &1.0315	&1.0343\\
         \hline
       	10   &1.1022	& 1.1037 &1.2302 &1.2365  & 1.03926	&1.04091 &1.23	&1.23\\
       	 \hline

	\varepsilonnd{tabular}
	\caption{Comparison between optimal Lagrange multipliers from Primal SDDP and Dual SDDP with penalties.}
	\label{tab-sen-num-dual}
\varepsilonnd{table}

With optimal dual solutions $\{\mbox{\boldmath$a$}r \pi_t\}$ and the realizations of $\{D_t\}$ and $\{{\varepsilonpsilon}silon_t\}$ at hand, we are able to compute the sensitivity of the optimal value with respect to $\phi$ and $\mu$, using (\ref{inv-12}) and (\ref{inv-13}),
with expectations estimated for $10\,000$ Monte Carlo simulations. We benchmark our method against the finite-difference  method. Specifically, for value function
$\vartheta$, the finite-difference method approximates the derivative with respect to $u_0$ by
$
	v'(u_0 ) \approx \frac{v(u_0 + \delta) - v(u_0 -\delta)}{2\delta}
	\label{sen-num:inv-3}$ for some small $\delta$.
	
	The sensitivity of the optimal value of the inventory problem with respect to $(\phi,\mu)$ is displayed in
	Table \ref{tab-sen-num-sensitivity}. In this table, S-$\phi$ and S-$\mu$ denote the derivatives with respect to
	$\phi$ and $\mu$ computed by our method, and fd-$\phi$, fd-$\mu$ denote the derivatives computed by the finite-difference method. In order to measure the difference between the two methods, we also compute S-gap-$\phi$ and S-gap-$\mu$, where
	S-gap-$\phi := \frac{|\text{fd-}\phi - \text{S-}\phi |}{|\text{fd-}\phi|} \times 100\%$ and S-gap-$\mu := \frac{|\text{fd-}\mu - \text{S-}\mu |}{|\text{fd-}\mu|} \times 100\%$.

 \begin{equation}gin{table}[htp]
 	\centering
 	\mathcal{S}mall
 	\begin{equation}gin{tabular}{|c|c|c|c|c|c|c|}
\hline
 		{\texttt{Instance}} & {\texttt{fd-}$\phi$} & {\texttt{S-}$\phi$} & {\texttt{S-gap-}$\phi$($\%$)}  & {\texttt{fd-}$\mu$} & {\texttt{S-}$\mu$} & {\texttt{S-gap-}$\mu$($\%$)}\\
 		\hline
 		1  & 403.604 & 401.094 & 0.622  & 164.578 & 164.158 & 0.255 \\
 		\hline
 		2   & 10\,716.111 & 10\,671.262 & 0.419 & 185.346 & 184.847 & 0.270 \\
 		\hline
 		3  & 269.514 & 269.443 & 0.026 & 134.646 & 134.463 & 0.136 \\
 		\hline
 		4   & 7\,780.570 & 7\,770.274 & 0.132 & 158.017 & 158.001 & 0.0101\\
 		\hline
 	\varepsilonnd{tabular}
 	\caption{Sensitivity of the optimal value with respect to $\phi$ and $\mu$ by the two methods.}
 	\label{tab-sen-num-sensitivity}
 \varepsilonnd{table}
We observe that the derivatives obtained by both methods are close to each other, especially when $\phi$ and $\mu$ are small. This is because small $\phi$ and $\mu$ gives rise to less variability in the demand. Note also that  the
finite-difference method is more time consuming since it requires computing the optimal value twice.
Instead, our method only needs to solve the model once. Moreover, computing the Lagrange multipliers does not significantly consume CPU time, as they are generated as a by-product of Primal SDDP. Alternatively, as discussed above, one can compute the optimal multipliers using Dual SDDP with penalties. Another drawback of the finite-difference method lies in its numerical instability. Indeed, the method is more accurate when $\delta$ is very small. However, the division by a very
small number generates bias while our approach is more stable.

\mathcal{S}ubsection{Dual SDDP for an hydro-thermal generation problem}

We repeat the experiments of Section \ref{sec:inventorysim} for the
Brazilian interconnected power system problem discussed  in \cite{pythonlibrarysddp} for $T=12$ stages
and $N_t=50$ inflow realizations for every stage. These realizations are
obtained calibrating log-normal distributions for each month of the year
using historical data of inflows and sampling from these distributions.
The data used for these simulations (including the inflow scenarios) is available on Github\footnote{{\tt{{\url{https://github.com/vguigues/Primal_SDDP_Library_Matlab}}}}}.

We solve this problem using Primal SDDP and Dual SDDP with penalizations. For this variant of Dual SDDP, a general
procedure to define sequences of penalizations
$(v_{t k})$ ensuring  convergence of the corresponding
Dual SDDP method is to take
$v_{t k}=\gamma_0 \alpha^{k-1}{\textbf{e}}$, $k \geq 1$, $t=2,\ldots,T$, with
$\alpha>1$, $\gamma_0>0$. For numerical reasons, we also take a large upper bound $U$ for these sequences and use
\begin{equation}gin{equation}\label{pensequence}
v_{t k} =  \min(U,\gamma_0 \alpha^{k-1}){\textbf{e}},\;k \geq 1, t=2,\ldots,T.
\varepsilonnd{equation}
We consider three variants of Dual SDDP: for the first variant, denoted by {\tt{Dual SDDP 1}},
$v_{t k}$  are as in \varepsilonqref{pensequence} with $\gamma_0=10^4$, $\alpha=1.3$, $U=10^{10}$.
To illustrate the fact that for constant sequences $v_{t k}=\gamma_0$, Dual SDDP converges (resp. does not converge)
for sufficiently large constants $\gamma_0$ (resp. sufficiently small constants $\gamma_0$)
we also define two other variants corresponding to  $U=+\infty$, $\gamma_0=10^9$, $\alpha=1$,
and $U=+\infty$, $\gamma_0=10^6$, $\alpha=1$, in \varepsilonqref{pensequence}, respectively denoted by
{\tt{Dual SDDP 2}} and {\tt{Dual SDDP 3}}.

We run Dual SDDP for 1000 iterations and Primal SDDP for 3000 iterations. The evolution of the upper and
lower bounds computed by the methods for the first 1000 iterations is given in Figure \ref{figboundsht}.\footnote{The upper bounds for Primal SDDP are computed as explained in Section \ref{sec:inventorysim}.}

\begin{equation}gin{figure}
 \centering
 \begin{equation}gin{tabular}{cc}
\includegraphics[scale=0.49]{BoundsT12M50HT.pdf}
 &
\includegraphics[scale=0.49]{BoundsT12M50HT1.pdf}\\
 \varepsilonnd{tabular}
\begin{equation}gin{tabular}{c}
\includegraphics[scale=0.5]{BoundsT12M50HT2.pdf}
 \varepsilonnd{tabular}
\caption{Top left: upper and lower bounds computed by Primal SDDP
and upper bounds computed by {\tt{Dual SDDP 1}}, {\tt{Dual SDDP 2}}, and {\tt{Dual SDDP 3}},
for the first 20 iterations for an instance of the hydro-thermal problem with $T=12$, $N_t=50$.
Top right: same outputs for iterations 21,$\ldots,150$.
Bottom: same outputs for iterations 151,$\ldots,1000$. }\label{figboundsht}
\varepsilonnd{figure}

More precisely, the values of these bounds for iterations 2, 5, 10, 50, 100, 150, 200, 250, 300, 350, 400, 1000,
and 3000 are reported in Table \ref{tablebounditerht}. We observe that parameter $\gamma_0$ for {\tt{Dual SDDP 3}}
is too small to allow this method to converge to the optimal value of the problem whereas the other two variants
{\tt{Dual SDDP 1}} and {\tt{Dual SDDP 2}} of Dual SDDP converge. Naturally, these methods start with
large upper bounds but after a few tens of iterations the upper bounds with {\tt{Dual SDDP 1}} and {\tt{Dual SDDP 2}}
are better than the upper bound computed by Primal SDDP. In particular, it is interesting to notice
that the best (lowest) upper bounds are obtained with the variant of Dual SDDP that uses adaptive
penalizations, i.e., penalizations that increase with the number of iterations before reaching value $U$
in \varepsilonqref{pensequence}.

\begin{equation}gin{table}
\centering
\begin{equation}gin{tabular}{|c|c|c|c|c|c|}
\hline
Iteration & \begin{equation}gin{tabular}{c}{\tt{Primal}}\\{\tt{SDDP}}\\Lb \varepsilonnd{tabular}&  \begin{equation}gin{tabular}{c}{\tt{Primal}}\\{\tt{SDDP}}\{\cal U}b\varepsilonnd{tabular}  & {\tt{Dual SDDP 1}} & {\tt{Dual SDDP 2}} & {\tt{Dual SDDP 3}}\\
\hline
2  & $1.317$ & $143.98$  & $1000.2$ & $1000.2$&$1000.2$\\
\hline
5  & $5.5588$ & $109.36$  & $1000.2$ &  $1000.2$&$994.04$\\
\hline
10 & $14.032$ &     $81.728$  & $360.40$   &  $1000.2$& $495.08$ \\
\hline
50     & $23.670$ &  $41.346$  & $54.999$  & $1000.2$& $96.720$ \\
\hline
100   & $24.787$ &  $35.502$  & $36.322$  & $64.072$ & $82.494$\\
\hline
150   & $25.111$ &  $32.447$  & $30.685$  & $35.595$& $79.465$\\
\hline
200   & $25.249$ &  $30.672$  & $29.076$  & $30.404$&$78.059$ \\
\hline
250   & $25.374$ &  $30.079$  & $28.215$  & $28.943$& $76.917$\\
\hline
300   & $25.436$ &  $29.434$  & $27.710$  & $28.030$& $76.344$\\
\hline
350   & $25.477$ &  $29.014$  & $27.309$  & $27.532$& $75.852$\\
\hline
400   & $25.526$ &  $28.626$  & $27.110$  & $27.188$& $75.526$\\
\hline
1000  & $25.703$ &  $27.175$  & $26.304$ & $26.335$  &$74.292$\\
\hline
3000  & $25.798$ &  $26.883$  & -  &  - &\\
\hline
\varepsilonnd{tabular}
\caption{For an instance of the hydro-thermal problem with $T=12$, $N_t=50$,
lower bound Lb and upper bound Ub computed by Primal SDDP and upper bounds computed
by variants of Dual SDDP along iterations. All costs have been divided
by $10^6$ to improve readability.}\label{tablebounditerht}
\varepsilonnd{table}

We also report in Table \ref{table1htiter} the relative error
$\frac{   {\tt{Upper}}_M(i)-{\tt{Lower}}_{{\tt{SDDP}}}(i)}{{\tt{Upper}}_M(i)}$ for iterations $i=100$, $200$, $300$, $400$, $500$, $800$, and $1000$
for all methods M where ${\tt{Upper}}_M(i)$ and ${\tt{Lower}}_{{\tt{SDDP}}}(i)$
are respectively the upper bound computed by method M at iteration $i$
and the lower bound computed by Primal SDDP at iteration $i$. For iterations 300 on, the relative error
is much smaller with variants of Dual SDDP, meaning that Primal SDDP overestimates the optimality gap.

\begin{equation}gin{table}
\centering
\begin{equation}gin{tabular}{|c|c|c|c|}
\hline
Iteration & {\tt{Primal SDDP}} & {\tt{Dual SDDP 1}} & {\tt{Dual SDDP 2}} \\
\hline
100  & 0.30 & 0.32 & 0.61\\
\hline
200  &  0.18& 0.13 &0.17\\
\hline
300  & 0.14  & 0.08 &0.09\\
\hline
400  &  0.11 & 0.06&0.06\\
\hline
500 &  0.09  & 0.05 &0.05\\
\hline
800 & 0.07  & 0.03 &0.03\\
\hline
1000 & 0.05   & 0.02 &0.02\\
\hline
\varepsilonnd{tabular}\\
\caption{Relative error as a function of the number of iterations for
{\tt{Primal SDDP}}, {\tt{Dual SDDP 1}}, and {\tt{Dual SDDP 2}}.
}\label{table1htiter}
\varepsilonnd{table}

However, each iteration of Dual SDDP takes more time as can be seen in
Figure \ref{fighttime} which reports the cumulative CPU time for all methods.
More precisely, running Dual and Primal SDDP in parallel, we can compute the time needed
to obtain a solution of relative accuracy ${\rm Var}epsilon$ using the standard stopping criterion
for Primal SDDP (see \cite{shapsddp}) or using the lower bound from Primal SDDP and the upper bound
from Dual SDDP, and computing the relative error obtained with these bounds each time a new bound
(either lower bound or upper bound) is computed. The results are reported in Table \ref{table1htprecision}.
We see that due to the fact that Dual SDDP iterations are more time consuming, for all relative accuracies but one,
the use of the stopping criterion based on Dual SDDP upper bounds requires more computational bulk.
From this experiment, performed on a larger problem (in terms
of size of the state vector and number of control variables for each stage) than the inventory problem
of Section \ref{sec:inventorysim}, it seems that the use of Dual SDDP
for a stopping criterion of Primal SDDP will decrease the overall computational
bulk only for small problems (having a limited to small number of controls, state variables, and
scenarios).

\begin{equation}gin{figure}
 \centering
\includegraphics[scale=0.65]{TimeHT.pdf}
\caption{Cumulative CPU time for Primal SDDP, Dual SDDP 1, Dual SDDP 2, and Dual SDDP 3.}\label{fighttime}
\varepsilonnd{figure}

\begin{equation}gin{table}
\centering
\begin{equation}gin{tabular}{|c|c|c|c|}
\hline
${\rm Var}epsilon$ & {\tt{Primal SDDP}} & {\tt{Dual SDDP 1}} & {\tt{Dual SDDP 2}} \\
\hline
0.3  & 515 & 1 042 &4 133\\
\hline
0.2  & 1 167 & 1 895 &7 446\\
\hline
0.15  & 1 659  & 2 910 &9 882\\
\hline
0.1  &  3 168 &5 114 &16 387\\
\hline
0.075 &  5 359  & 8 003 &22 457\\
\hline
0.05 &  11 124  & 15 738 &35 113\\
\hline
0.04 &  45 391  & 23 449  &51 381\\
\hline
\varepsilonnd{tabular}\\
\caption{Time (in seconds) needed to obtain a solution of relative accuracy ${\rm Var}epsilon$
with Primal SDDP and variants of Dual SDDP for an instance of the hydro-thermal problem.}\label{table1htprecision}
\varepsilonnd{table}

\if{
Finally, similarly to Figure \ref{fig3error},
as an evidence of the fact that RCR still does not hold for the dual of the hydro-thermal problem, we
observed that the maximal
 and mean values
 of  $\|{ z}eta_t^{k}\|_1$
along iterations, where ${ z}eta_t^{k}$ is an optimal
value
of ${ z}eta_t$ in \varepsilonqref{forwarddualsddppen} for iteration $k$, are positive for some stages.
}\fi

\addcontentsline{toc}{section}{References}

\mathcal{S}etcounter{equation}{0}
\mathcal{S}ection{Appendix}
\label{appen}
In this Appendix, we prove
Lemma \ref{boundedpi}, Proposition \ref{proofubfeassddp}, and Theorem \ref{thconvdualsddppenal}.

We first need more notation.
We introduce the sequence of  functions
\begin{equation}gin{equation}\label{defvhatTbis}
{\overline V}_T^k (\pi_{T-1}):=
\begin{equation}gin{array}{l}
\max\limits_{\pi_{T 1},\ldots,\pi_{T N_T}, { z}eta_T} \; \displaystyle \mathcal{S}um\limits_{j=1}^{N_T} p_{T j}
 b_{T j}^\top \pi_{T j}-v_{T k}^\top { z}eta_T\\
A_{T}^{\top} \pi_{T j} \leq  c_{T},\;j=1,\ldots,N_T,\\
A_{T-1}^{\top} \pi_{T-1} + \displaystyle \mathcal{S}um_{j=1}^{N_T}
p_{T j} B_{T j}^{\top} \pi_{T j} \leq c_{T-1}
+ { z}eta_T,\\
{ z}eta_T \geq 0, {\underline \pi}_T \leq \pi_{T j} \leq {\overline \pi}_{T},\;j=1,\ldots,N_T,
\varepsilonnd{array}
\varepsilonnd{equation}
and for $t=2,\ldots,T-1$, the sequence of functions
\begin{equation}gin{equation}\label{defvhattbis}
{\overline V}_t^k(\pi_{t-1}):=
\begin{equation}gin{array}{l}
\max\limits_{\pi_{t 1},\ldots,\pi_{t N_t}, { z}eta_t} \; \displaystyle \mathcal{S}um\limits_{j=1}^{N_t} p_{t j} \left(
  b_{t j}^\top \pi_{t j} + {V}_{t+1}^k (\pi_{t j}) \right) - v_{t k}^\top { z}eta_t\\
A_{t-1}^{\top} \pi_{t-1}+ \displaystyle \mathcal{S}um_{j=1}^{N_t} p_{t j} B_{t j}^{\top} \pi_{t j}  \leq  c_{t-1} + { z}eta_t,\\
{ z}eta_{t} \geq 0,{\underline \pi}_t \leq \pi_{t j} \leq {\overline \pi}_{t},\;j=1,...,N_t.\\
\varepsilonnd{array}
\varepsilonnd{equation}
Due to Assumption (A1) we can
represent the scenarios
for ${ x}i_1,{ x}i_2,\ldots,{ x}i_T$,
by a scenario tree
of depth $T+1$
where the root node $n_0$ associated to a stage $0$ (with decision $x_0$ taken at that
node) has one child node $n_1$
associated to the first stage. We denote by 
$\mathcal{N}$ the set of nodes and for a node $n$ 
of the tree, by $(x_n,\pi_n)$ a primal-dual pair at that 
node and by
${ x}i_n$ the realization of process $({ x}i_t)$ at node $n$ (this realization ${ x}i_n$ contains in particular the realizations
$c_n$ of $c_t$, $b_n$ of $b_t$, $A_{n}$ of $A_{t}$, and $B_{n}$ of $B_{t}$).\\

\par {\textbf{Proof of Lemma \ref{boundedpi}.}}
Let $1 \leq t \leq T$ and
let us fix a node $m$
of stage $t$.
Let
$\overline A_m$
such that
constraints
$A_m x_m + B_m x_{F(m)}=b_m$
are rewritten in compact form
$\overline A_m x = b_m$
in terms of vector
$x=(x_n)_{n \in \mathcal{N}}$
of decision variables
in the scenario tree.
The dual function obtained dualizing
the coupling constraints of node
$m$ is given by
$$
\theta(\pi_m) =
\begin{equation}gin{array}{l}
\min \; \mathbb{E}[c^\top x] + \pi_m^\top
(A_m x_m + B_m x_{F(m)}-b_m)\\
x \in \mathcal{S}_m,
\varepsilonnd{array}
$$
for
$
\mathcal{S}_m
=\{
x=(x_n)_{n \in \mathcal{N}} :
x \geq 0 \} \cap \mathcal{A}_m
$
where
$\mathcal{A}_m=\{
x=(x_n)_{n \in \mathcal{N}} :
 A_n x_n +B_n x_{F(n)} = b_n, \forall
n \neq m, n \in \mathcal{N} \}.
$

By Linear Programming Duality, the optimal
value $\mathcal{Q}_1(x_0)$ of primal
problem \varepsilonqref{eq-1}
is the optimal value of the dual
problem
\begin{equation}gin{equation}\label{dualpbinit}
\max \{\theta(\pi_m) : \pi_m \in \mathbb{R}^{m_t}\}
\varepsilonnd{equation}
which can clearly be written as
\begin{equation}gin{equation}\label{reformdual}
\mathcal{Q}_1(x_0) = \max_{\pi_m}\{ \theta(\pi_m):
\pi_m = {\overline A}_m x -
{b}_m, x \in \mbox{Aff}(\mathcal{S}_m)\},
\varepsilonnd{equation}
where $\mbox{Aff}(\mathcal{S}_m)$
is the affine hull of  $\mathcal{S}_m$.
We now bound the optimal solutions
of dual problem \varepsilonqref{reformdual}.
Since \varepsilonqref{dualpbinit}
and \varepsilonqref{reformdual} have the same
optimal values, adding these bounds as
constraints on $\pi_m$ in
\varepsilonqref{dualpbinit} does not change
its optimal value. Since $\hat x>0$ there is $r>0$ such that
\begin{equation}gin{equation}\label{xhatrelint}
\mathbb{B}(\hat x,r)
\mathcal{S}ubseteq \{x \geq 0\}.
\varepsilonnd{equation}
We argue that
$\mbox{Aff}(\mathcal{S}_m)=\mathcal{A}_m$.
Indeed, the inclusion
$\mbox{Aff}(\mathcal{S}_m) \mathcal{S}ubseteq \mathcal{A}_m$ is clear.
Now if $x \in \mathcal{A}_m$ then if $x = \hat x$
we have that $x \in \mathcal{S}_m \mathcal{S}ubseteq
\mbox{Aff}(\mathcal{S}_m)$ and if $x \neq \hat x$,
recalling
that $\hat x \in \mathcal{A}_m$ satisfies
\varepsilonqref{xhatrelint} we have that
$$
y:=\hat  x + \frac{r}{2} \frac{x-\hat x}{\|x-\hat x\|}
\in \mathcal{A}_m \cap \mathbb{B}(\hat x,r)
\mathcal{S}ubseteq \mathcal{S}_m.
$$
Therefore $x$ belongs to the line that contains $y$ and $\hat x$ with
$y, \hat x$ belonging to
$\mathcal{S}_m$ which implies $x \in \mbox{Aff}(\mathcal{S}_m)$
and
$\mbox{Aff}(\mathcal{S}_m)=\mathcal{A}_m$.

It follows that
$$
\mathbb{B}(\hat x,r)
\cap \mbox{Aff}(\mathcal{S}_m)
= \mathbb{B}(\hat x,r)
\cap \mathcal{A}_m \mathcal{S}ubseteq \mathcal{S}_m
$$
and that
 there is $\rho_*(m)>0$
such that
$$
\mathbb{B}(0,\rho_*)
\cap (\overline{A}_m \mathcal{A}_m - b_m)
\mathcal{S}ubseteq \overline A_m
(\mathbb{B}(\hat x,r) \cap \mathcal{A}_m)
-b_m.
$$
Let $\mbox{\boldmath$a$}r \pi_m$ be
an optimal solution of problem \varepsilonqref{reformdual}
and let $z=0$ if $\mbox{\boldmath$a$}r \pi_m=0$
and
$z=-\frac{\mbox{\boldmath$a$}r \pi_m}{\|\mbox{\boldmath$a$}r \pi_m\|_2}\rho_*$
otherwise.
Observe that
$z \in \mathbb{B}(0,\rho_*)
\cap (\overline{A}_m \mathcal{A}_m - b_m)
$
and therefore
$z \in \overline A_m
(\mathbb{B}(\hat x,r) \cap \mathcal{A}_m)
-b_m \mathcal{S}ubseteq
\overline A_m \mathcal{S}_m-b_m$
and $z$ can be written
$z={\overline A}_m \tilde x - b_m$
for $\tilde x \in \mathbb{B}(\hat x,r) \cap \mathcal{S}_m$. It follows
that
$$
\begin{equation}gin{array}{lcl}
\mathcal{Q}_1(x_0)
=\theta(\mbox{\boldmath$a$}r \pi_m)
& \leq &
\mathbb{E}[c^\top \tilde x] +
\mbox{\boldmath$a$}r \pi_m^{\top}({\overline A}_m \tilde x - b_m)\\
& \leq & \mathbb{E}[c^\top \hat x] + r
\mathcal{S}um_{t=1}^T \mathbb{E}[\| c_t \|_2] +  \mbox{\boldmath$a$}r \pi_m^{\top} z\\
& = &  \mathbb{E}[c^\top \hat x] +
r \mathcal{S}um_{t=1}^T \mathbb{E}[\| c_t \|_2]  - \rho_*(m) \|\mbox{\boldmath$a$}r \pi_m\|_2
\varepsilonnd{array}
$$
which gives for every node
$n$ of stage $t$ that
$$
\|\mbox{\boldmath$a$}r \pi_n\|_2
\leq
\max_{m \in {\tt{Nodes}}(t)}
\frac{
\mathbb{E}[c^\top \hat x]
-\mathcal{Q}_1(x_0)+r \mathcal{S}um_{t=1}^T \mathbb{E}[\| c_t \|_2]}{\rho_*(m)}
$$
with corresponding box constraints
$\underline \pi_t, \overline \pi_t$
where ${\tt{Nodes}}(t)$ are the nodes
of stage $t$. $
 \mathcal{S}quare$\\

\par {\textbf{Proof of Proposition \ref{proofubfeassddp}.}} We show by induction on $k$ that $V_t \leq {V}_t^k$ for $t=2,\ldots,T$. For $k=0$ these relations hold by definition. Assume that for some $k \geq 1$ we have $V_t \leq {V}_t^{k-1}$ for $t=2,\ldots,T$. We show by backward induction on $t$ that $V_t \leq {V}_t^{k}$ for $t=2,\ldots,T$. Observe that for any $\pi_{T-1}$, optimization problem \varepsilonqref{defvhatTbis} with optimal value ${\overline V}_T^k(\pi_{T-1})$ is feasible. Indeed, since primal problem \varepsilonqref{eq-1} is feasible and has a finite optimal value, the corresponding dual problem is feasible which implies that there is $\pi_{T 1},\ldots,\pi_{T N_T}$ satisfying $A_{T}^{\top} \pi_{T j} \leq  c_{T}$, ${\underline \pi}_T \leq \pi_{T j} \leq {\overline \pi}_{T}$, $j=1,\ldots,N_T$, and for every such points we can find ${ z}eta_T \geq 0$ satisfying the remaining constraints in \varepsilonqref{defvhatTbis}. Therefore ${\overline V}_T^k(\pi_{T-1})$ is finite for every $\pi_{T-1}$ and is the optimal value of the corresponding dual optimization problem, i.e., for any $\pi_{T-1}$ we get $$ {\overline V}_T^k(\pi_{T-1}) = \begin{equation}gin{array}{l} \displaystyle \min_{\alpha, \delta, {\overline{\Psi}}, {\underline{\Psi}}} \; \delta^{\top} (c_{T-1} - A_{T-1}^{\top} \pi_{T-1})  + c_T^{\top}  \mathcal{S}um_{j=1}^{N_T} \alpha_j + \mathcal{S}um_{j=1}^{N_T} {\overline{\Psi}}_j^\top {\overline \pi}_T - \mathcal{S}um_{j=1}^{N_T} {\underline{\Psi}}_j^\top {\underline \pi}_T \\ A_T \alpha_j + p_{T j} B_{T j} \delta - {\underline{\Psi}}_j +{\overline{\Psi}}_j = p_{T j} b_{T j},\;j=1,\ldots,N_T,\\ 0 \leq \delta \leq v_{T k}, \alpha_j, {\underline{\Psi}}_j, {\overline{\Psi}}_j \geq 0,\;j=1,\ldots,N_T. \varepsilonnd{array} $$ Using this dual representation and the definition of $\overline{\theta}_T^k, {\overline{\begin{equation}ta}}_T^k$, we get for every $\pi_{T-1}$: \begin{equation}gin{equation}\label{cutaboveT} \overline{\theta}_T^k + \langlegle {\overline{\begin{equation}ta}}_T^k , \pi_{T-1} {\ranglegle}ngle  \geq {\overline V}_T^k(\pi_{T-1}). \varepsilonnd{equation} Recalling representation \varepsilonqref{defvhatTbis} for ${\overline V}_T^k(\pi_{T-1})$, observe that for every $\pi_{T-1} \in \mbox{dom}(V_T)$ we have ${\overline V}_T^k(\pi_{T-1}) \geq V_T(\pi_{T-1})$ whereas for $\pi_{T-1} \notin \mbox{dom}(V_T)$ we have $V_T(\pi_{T-1})=-\infty$ while ${\overline V}_T^k(\pi_{T-1})$ is finite, which shows that for every $\pi_{T-1}$ we have $ {\overline V}_T^k(\pi_{T-1}) \geq V_T(\pi_{T-1})$, which, combined with \varepsilonqref{cutaboveT} and the induction hypothesis, gives $$V_T^k(\pi_{T-1}) \geq V_T(\pi_{T-1})$$ for every$\pi_{T-1}$.

Now assume that $V_{t+1}^k(\pi_{t}) \geq V_{t+1}(\pi_{t})$ for all $\pi_t$ for some  $t \in \{2,\ldots,T-1\}$. We want to show that $V_{t}^k(\pi_{t-1}) \geq V_{t}(\pi_{t-1})$ for all $\pi_{t-1}$. First observe that for every $\pi_{t-1}$, linear program \varepsilonqref{defvhattbis} with optimal value ${\overline V}_t^k(\pi_{t-1})$ is feasible and has a finite optimal value. Therefore we can express ${\overline V}_t^k(\pi_{t-1})$ as the optimal value of the corresponding dual problem given by \begin{equation}gin{equation}\label{dualpbvtk1} \begin{equation}gin{array}{l} \displaystyle \min_{\delta,\nu,{\overline \Psi}, {\underline \Psi}} \; \delta^{\top}{\cal B}ig[ c_{t-1} - A_{t-1}^{\top} \pi_{t-1}{\cal B}ig] + \mathcal{S}um_{i=0}^{k} {\overline{\theta}}_{t+1}^i \mathcal{S}um_{j=1}^{N_t} \nu_i(j) + \mathcal{S}um_{j=1}^{N_t} {\overline \Psi}_j^\top {\overline \pi}_t -  \mathcal{S}um_{j=1}^{N_t} {\underline \Psi}_j^\top {\underline \pi}_t\\ \displaystyle p_{t j} B_{t j} \delta - \mathcal{S}um_{i=0}^k \nu_i(j)  {\overline{\begin{equation}ta}}_{t+1}^i - {\underline \Psi}_j + {\overline \Psi}_j = p_{t j} b_{t j},\;j=1,\ldots,N_t,\\ \displaystyle \mathcal{S}um_{i=0}^k \nu_i( j ) = p_{t j}, {\underline \Psi}_j, {\overline \Psi}_j \geq 0,\;j=1,\ldots,N_t,\\ \nu_0,\ldots,\nu_k \geq 0,0 \leq \delta \leq v_{t k}. \varepsilonnd{array} \varepsilonnd{equation} Using this representation of $\overline V_t^k$ and the definition of $\overline{\theta}_t^k, {\overline{\begin{equation}ta}}_t^k$, we obtain for every $\pi_{t-1}$: \begin{equation}gin{equation}\label{cutabovet} \overline{\theta}_t^k +  \langlegle {\overline{\begin{equation}ta}}_t^k , \pi_{t-1} {\ranglegle}ngle \geq {\overline V}_t^k(\pi_{t-1}). \varepsilonnd{equation} Next, recalling representation \varepsilonqref{defvhattbis} for ${\overline V}_t^k(\pi_{t-1})$ and the induction hypothesis, we get \begin{equation}gin{equation}\label{induction2} {\overline V}_t^k(\pi_{t-1}) \geq {\omegaidehat V}_t^k(\pi_{t-1}) \varepsilonnd{equation} where $$ {\omegaidehat V}_t^k(\pi_{t-1}):= \begin{equation}gin{array}{l} \max\limits_{\pi_{t 1},\ldots,\pi_{t N_t}, { z}eta_t} \; \displaystyle \mathcal{S}um\limits_{j=1}^{N_t} p_{t j} \left( b_{t j}^\top \pi_{t j} + {V}_{t+1}(\pi_{t j}) \right) - v_{t k}^\top { z}eta_t\\ A_{t-1}^{\top} \pi_{t-1}+ \displaystyle \mathcal{S}um_{j=1}^{N_t} p_{t j} B_{t j}^{\top} \pi_{t j}  \leq  c_{t-1} + { z}eta_t,\\ { z}eta_{t} \geq 0,{\underline \pi}_t \leq \pi_{t j} \leq {\overline \pi}_{t},\;j=1,...,N_t.\\ \varepsilonnd{array} $$ Similarly to the induction step $t=T$, for every $\pi_{t-1}$, we have \begin{equation}gin{equation}\label{induction3} {\omegaidehat V}_t^k(\pi_{t-1}) \geq V_t(\pi_{t-1}). \varepsilonnd{equation} Combining \varepsilonqref{cutabovet}, \varepsilonqref{induction2}, and \varepsilonqref{induction3} with the induction hypothesis, we obtain $V_{t}^k(\pi_{t-1}) \geq V_{t}(\pi_{t-1})$ for all $\pi_{t-1}$ which achieves the proof of the induction step $t$.

In particular $V_2^{k-1} \geq V_2$ which implies that $V^{k-1}$ is greater than or equal to the optimal value of dual problem \varepsilonqref{eq-2} which is also, by linear programming duality, the optimal value of primal problem \varepsilonqref{eq-1}.$
 \mathcal{S}quare$\\

The proof of Theorem \ref{thconvdualsddppenal} is based on the following lemma:

\begin{equation}gin{lemma}
\label{lem-conver}
Suppose that the multistage  problem  \varepsilonqref{eq-1} has a finite optimal value. Then for sufficiently large values of the components of vectors $v_t$,  in the dynamic equations \varepsilonqref{eq-r2}, the optimal   value of the multistage problem defined by these dynamic equations coincides with the optimal value of the original problem \varepsilonqref{eq-1}.
\varepsilonnd{lemma}

 {\bf Proof.} 
  As it was already mentioned,  since it is assumed that the number of scenarios is finite,
 we can view problem \varepsilonqref{eq-1} as a large  linear program (deterministic equivalent) written under the form 
 \begin{equation}gin{equation}\label{linprimalpr-1}
\min_{x} c^\top x\;\;{\rm s.t.}\; 
{\cal A} x =b,\;x \geq 0.  
 \varepsilonnd{equation}
Also since \varepsilonqref{eq-1} has a finite optimal
 value, it has a nonempty set of optimal solutions
 and there is a bounded optimal solution 
of \varepsilonqref{linprimalpr-1}.
Let us fix such an optimal solution $\mbox{\boldmath$a$}r x$. 
We have that problem \varepsilonqref{linprimalpr-1}  can
 be written
 \begin{equation}gin{equation}\label{linprimalpr-2}
\min_{x} c^\top x\;\;{\rm s.t.}\; 
{\cal A} x =b,\;0 \leq x \leq \mbox{\boldmath$a$}r x.  
 \varepsilonnd{equation}
 The dynamic programming equations \varepsilonqref{eq-3a-1} - \varepsilonqref{eq-7-1} represent the standard dual of \varepsilonqref{eq-1}. We can also  think about that dual as a large linear programming problem of the form
 (this is the dual of \varepsilonqref{linprimalpr-1}):
 \begin{equation}gin{equation}\label{lindualpr-1}
\max_{\pi} b^\top \pi\;\;{\rm s.t.}\; 
{\cal A}^\top \pi \le c.  
 \varepsilonnd{equation}
Similarly the deterministic equivalent of  penalized dynamic equations \varepsilonqref{eq-r2} can be written as:
 \begin{equation}gin{equation}\label{lindualpr-2}
\max_{\pi,{ z}eta} b^\top  \pi-v^\top { z}eta\;\;{\rm s.t.}\;
{\cal A}^\top \pi \le c+{ z}eta, { z}eta \geq 0.
 \varepsilonnd{equation}
Next, from optimality conditions of linear programs,
$(x,\pi)$ is an optimal primal-dual pair
for \varepsilonqref{linprimalpr-1}-\varepsilonqref{lindualpr-1}
if and only if 
 \begin{equation}gin{equation}\label{lindualpr-3}
x^\top({\cal A}^\top \pi -c)=0,\;
{\cal A} x =b,\;x \geq 0,\; 
{\cal A}^\top \pi \leq c.
 \varepsilonnd{equation}
 The corresponding optimality conditions for \varepsilonqref{lindualpr-2} are
 \begin{equation}gin{equation}\label{lindualpr-4}
 x^\top({\cal A}^\top \pi-c-{ z}eta)-{ z}eta^\top \gamma=0,\;
{\cal A}^\top \pi \leq c+ { z}eta,\;
{ z}eta \geq 0,\;
{\cal A} x = b,\;x \geq 0,\;\gamma \geq 0,\;
x=v-\gamma. 
\varepsilonnd{equation}
Now let $\mbox{\boldmath$a$}r{\pi}$ be an optimal dual solution, i.e.,
an optimal solution of \varepsilonqref{lindualpr-1}.
Then \varepsilonqref{lindualpr-3}
is satisfied with 
$(x,\pi)=(\mbox{\boldmath$a$}r x, \mbox{\boldmath$a$}r \pi)$.
It follows  that if $v\ge \mbox{\boldmath$a$}r x$, then 
$(x,\pi,{ z}eta,\gamma)=(\mbox{\boldmath$a$}r x,\mbox{\boldmath$a$}r \pi,0,v-\mbox{\boldmath$a$}r x)$  with  
${ z}eta=0$ satisfies \varepsilonqref{lindualpr-4}, and hence 
$(\mbox{\boldmath$a$}r \pi,\mbox{\boldmath$a$}r { z}eta)=(\mbox{\boldmath$a$}r \pi,0)$
is an optimal solution of  
\varepsilonqref{lindualpr-2} showing that the optimal value
of \varepsilonqref{lindualpr-2} is
$b^\top \mbox{\boldmath$a$}r \pi=c^\top \mbox{\boldmath$a$}r x$, i.e., the
optimal value of \varepsilonqref{linprimalpr-1}.
We obtain  that 
for $v\ge \mbox{\boldmath$a$}r x$, the optimal values of problems  \varepsilonqref{lindualpr-1} and \varepsilonqref{lindualpr-2} do coincide.\footnote{Observe that the dual of \varepsilonqref{lindualpr-2}
is given by
$$
\min_{x} c^\top x\;\;{\rm s.t.}\; 
{\cal A} x =b,\;0 \leq x \leq v,
$$
and for $v \geq \mbox{\boldmath$a$}r x$, this linear program has the same
optimal value as \varepsilonqref{linprimalpr-2}, which, as we have seen, is equivalent to primal problem \varepsilonqref{eq-1}.}$
 \mathcal{S}quare$\\

\par {\textbf{Proof of Theorem \ref{thconvdualsddppenal}.}} Dual SDDP with
penalizations is SDDP applied to Dynamic Programming equations  
corresponding to a linear program with finite optimal value, satisfying relatively complete recourse
with discrete uncertainties of finite support. Since samples $\tilde { x}i_t^k$
in Dual SDDP with penalizations are independent, 
we can follow the convergence proof of SDDP for linear programs
from \cite{philpot} 
to obtain that  $V^k$ converges
to the optimal value of the penalized linear programs, which,
by Lemma \ref{lem-conver} (observe that the Lemma can be applied
since $\lim_{k \rightarrow +\infty} v_{t k}=+\infty$),
is the optimal value  of \varepsilonqref{eq-1}.$
 \mathcal{S}quare$

\if{

Let us take an arbitrary realization of the algorithm and let us show that
$V^k$ converges to
the optimal value of
\varepsilonqref{eq-1}
for that realization. Define $V_{T+1}$ and
${\overline{V}}_{T+1}^k$ for $k \geq 0$
on $[\underline{\pi}_T,\overline{\pi}_T]$
to be the null function.
We show that
\begin{equation}gin{equation}\label{inductionkforward}
\lim_{k \rightarrow +\infty} V_t^k(\pi_n^k)-V_t(\pi_n^k)=0
\varepsilonnd{equation}
almost surely
for all $t=2,\ldots,T+1$ and all node $n$ of stage
$t-1$ where $\pi_n^k$ is computed simulating the dual policy obtained
using Dual SDDP with penalizations replacing $V_t$ by $V_t^{k-1}$.
The conclusion of the theorem is an immediate consequence of
\varepsilonqref{inductionkforward} written for $t=2$.

We show \varepsilonqref{inductionkforward} by backward induction on $t$.
The relation holds for $t=T+1$. Assume it holds for
$t+1$ with $t \in \{2,\ldots,T\}$ and let us show it holds for stage $t$.
Take a node $n$ of stage $t-1$.
It is enough to show that
$\lim_{k \rightarrow +\infty, k \in \mathcal{S}_n} V_t^k(\pi_n^k)-V_t(\pi_n^k)=0$
where $\mathcal{S}_n=\{k: ({ x}i_1,\tilde { x}i_2^k,\ldots,\tilde { x}i_{t-1}^k) =
{ x}i_{[n]}\}$ with ${ x}i_{[n]}$ the scenario from node $n_1$
to node $n$. Indeed, as in \cite{lecphilgirar12, guiguessiopt2016}, from this relation, \varepsilonqref{inductionkforward}
can then be shown by contradiction using the
Strong Law of Large Numbers and the assumption of independence
of ${ x}i_t^k$. Therefore let $k \in \mathcal{S}_n$.
Since $V_t^{k} \leq V_t^{k-1}$ (by definition
of $V_t^k)$
and $v_{t k}\geq v_{t k-1}$,
we have
\begin{equation}gin{equation}\label{decVtkbar}
{\overline V}_t^k \leq {\overline V}_t^{k-1}.
\varepsilonnd{equation}
Since $k \in \mathcal{S}_n$, we have
$\overline \theta_t^k
+ \langlegle \overline \begin{equation}ta_t^k , \pi_{n}^k {\ranglegle}ngle
= \overline V_t^k(\pi_{n}^k)$ and therefore
for $0 \leq j \leq k$ we get
$
\overline \theta_t^k
+ \langlegle \overline \begin{equation}ta_t^k , \pi_{n}^k {\ranglegle}ngle
= \overline V_t^k(\pi_{n}^k)
 \leq  {\overline V}_t^j(\pi_{n}^k)
\leq \overline \theta_t^j
+ \langlegle \overline \begin{equation}ta_t^j , \pi_{n}^k {\ranglegle}ngle,
$
which implies, by definition of $V_t^k$,
\begin{equation}gin{equation}\label{decVtkbar2}
V_t^k(\pi_{n}^k) = {\overline V}_t^k(\pi_{n}^k).
\varepsilonnd{equation}
Therefore for $k \in \mathcal{S}_n$, we have
\begin{equation}gin{equation}\label{proofconvdsddp}
0  \leq  V_t^k(\pi_n^k)-V_t(\pi_n^k)    \mathcal{S}tackrel{\varepsilonqref{decVtkbar2}}{=}  {\overline V}_t^k(\pi_n^k)-V_t(\pi_n^k)
 \leq  {\overline V}_t^{F, k-1}(\pi_n^k)-V_t(\pi_n^k)
\varepsilonnd{equation}
where ${\overline V}_t^{F, k-1}(\pi_n^k)$ is the optimal value of
\begin{equation}gin{equation}\label{VFtk}
\begin{equation}gin{array}{l}
\max\limits_{\pi_{m}, { z}eta_t} \; \displaystyle \overline F_t^{k-1}((\pi_m)_{m \in C(n)})= \mathcal{S}um\limits_{m \in C(n)} p_{m} \left(
  b_{m}^\top \pi_{m} + {V}_{t+1}^{k-1} (\pi_{m}) \right) - v_{t k}^\top { z}eta_t\\
A_{t-1}^{\top} \pi_{n}^k + \displaystyle \mathcal{S}um_{m \in C(n)} p_{m} B_{m}^{\top} \pi_{m}  \leq  c_{t-1} + { z}eta_t,\\
{ z}eta_{t} \geq 0,{\underline \pi}_t \leq \pi_{m} \leq {\overline \pi}_{t},\;m \in C(n).\\
\varepsilonnd{array}
\varepsilonnd{equation}
The problem above is the problem
solved in the forward pass for
iteration $k$ stage $t$.
By first order optimality conditions,
$(\pi^k=(\pi_m^k)_{m \in C(n)},{ z}eta_t^k)$ is
an optimal primal solution of this problem
and $(\delta_t^k, \lambda_t^k, {\underline \Psi}_m^k, {\overline \Psi}_m^k)$\footnote{We use for dual variables the same notation as in dual problem \varepsilonqref{dualpbvtk1}} is an optimal dual
solution if and only they solve the following
KKT system:
\begin{equation}gin{equation}\label{kkt1}
\begin{equation}gin{array}{l}
A_{t-1}^\top \pi_n^k + \mathcal{S}um_{m \in C(n)} p_m
B_n^{\top} \pi_m^k \leq c_{t-1}+{ z}eta_t^k, { z}eta_t^k \geq 0, \underline \pi_t^k \leq \pi_m^k \leq \overline \pi_t,\forall m \in C(n),\\
\delta_t^k, \lambda_t^k, \underline \Psi_m^k, \overline \Psi_m^k \geq 0,\\
0=-p_m(b_m+s_m)+p_mB_m \delta_t^k - \underline \Psi_m^k
+ \overline \Psi_m^k,\\
0=v_{t k} -\delta_t^k - \lambda_t^k,\\
\lambda_t^k(i){ z}eta_t^k(i)=0,\underline{\Psi}_m(i)
(\underline{\pi}_t(i)-\pi_m(i))=0,
\overline{\Psi}_m(i) (\pi_m(i)-\overline{\pi}_t(i))=0,\forall i,\\
\delta_t^k(i)(A_{t-1}^\top \pi_n^k + \mathcal{S}um_{m \in C(n)} p_m
B_n^{\top} \pi_m^k - c_{t-1}-{ z}eta_t^k)(i)=0, \forall i,
\varepsilonnd{array}
\varepsilonnd{equation}
for some $s_m \in \partial V_{t+1}^{k-1}(\pi_m^k)$.
Let
$$\mathcal{X}_t(\pi_n)=
\{(\pi_m):
A_{t-1}^{\top} \pi_{n} + \mathcal{S}um_{m \in C(n)} p_{m} B_{m}^{\top} \pi_{m}  \leq  c_{t-1}, {\underline \pi}_t \leq \pi_{m} \leq {\overline \pi}_{t},\;m \in C(n)\}.$$
Let also $({\hat \pi}_m^k)_{m \in C(n)}$
be an optimal solution of
\begin{equation}gin{equation}\label{VFtkbis}
\begin{equation}gin{array}{l}
\max\limits_{\pi_{m}} \; \displaystyle
F_t^{k-1}((\pi_m)_{m \in C(n)})=
 \mathcal{S}um\limits_{m \in C(n)} p_{m} \left(
  b_{m}^\top \pi_{m} + {V}_{t+1}^{k-1} (\pi_{m}) \right)\\
(\pi_m)_{m \in C(n)} \in  \mathcal{X}_t(\pi_n)
\varepsilonnd{array}
\varepsilonnd{equation}
and let
$(\hat \delta_t^k, {\hat{{\underline \Psi}}}_m^k,
{\hat{{\overline \Psi}}}_m^k)$
be an optimal dual solution (using a notation similar as before
for dual variables).
We partition $\mathcal{S}_n$ in $D_n \cup {\overline{D}}_n$
where $D_n$ is the set of iterations in
$\mathcal{S}_n$ such that
$F_t(\pi_n^k)$ is nonempty.
We first argue that due to
the compactness of
$[\underline{\pi}_t,\overline{\pi}_t]$ and the equicontinuity of
$V_t^k$ there can only be a finite
number of iterations in $\overline D_n$
so that after some  iteration $k_0$
all iterations in $\mathcal{S}_n$
are in $D_n$.
Next, following the proof
of Lemma \ref{boundedpi}, using the fact that
$\underline{\pi}_t<\overline{\pi}_t$,
and that $F_t^{k-1}\leq F_t^{0}$
with $F_t^{0}$ Lipschitz continuous, we get that
there is $\delta_*$  such that for all
$k$ we have $0 \leq \hat \delta_t^k \leq \delta_*$.
From the first order optimality conditions,
$({\hat \pi}_m^k)_{m \in C(n)}, (\hat \delta_t^k, {\hat{{\underline \Psi}}}_m^k,
{\hat{{\overline \Psi}}}_m^k)$
is an optimal primal-dual solution
of \varepsilonqref{VFtkbis} if and only if
they solve the following
KKT system:
\begin{equation}gin{equation}\label{kkt2}
\begin{equation}gin{array}{l}
A_{t-1}^\top \pi_n^k + \mathcal{S}um_{m \in C(n)} p_m
B_n^{\top} \hat \pi_m^k \leq c_{t-1}, \underline \pi_t^k \leq \hat \pi_m^k \leq \overline \pi_t,\forall m \in C(n),\\
\hat \delta_t^k, {\hat{\underline \Psi}}_m^k,
{\hat{\overline{\Psi}}}_m^k \geq 0,\\
0=-p_m(b_m+s_m)+p_m B_m {\hat \delta}_t^k - {\hat{\underline \Psi}}_m^k
+ {\hat{\overline \Psi}}_m^k,\\
(\underline{\pi}_t(i)-\hat \pi_m(i))=0,
{\hat{\overline{\Psi}}}_m(i) ({\hat{\pi}}_m(i)-\overline{\pi}_t(i))=0,\forall i,\\
{\hat{\delta}}_t^k(i)(A_{t-1}^\top \pi_n^k + \mathcal{S}um_{m \in C(n)} p_m
B_n^{\top} \hat \pi_m^k - c_{t-1})(i)=0, \forall i,
\varepsilonnd{array}
\varepsilonnd{equation}
for some $s_m \in \partial V_{t+1}^{k-1}(\pi_m^k)$.
Now due to $\lim_{k \rightarrow +\infty} v_{tk}=+\infty$
we can choose $k_0$ such that
for $k \geq k_0$ we have
$v_{t k}>\delta_*$ and therefore
any primal-dual optimal solution $(
({\hat \pi}_m^k)_{m \in C(n)}, (\hat \delta_t^k, {\hat{{\underline \Psi}}}_m^k,
{\hat{{\overline \Psi}}}_m^k))$
of \varepsilonqref{VFtkbis} which satisfies
\varepsilonqref{kkt2} provides a primal-dual optimal solution
for \varepsilonqref{VFtk} given by
$(\pi_m^k)_{m \in C(n)}=({\hat \pi}_m^k)_{m \in C(n)}$,
${ z}eta_t^k=0$, $\lambda_t^k=v_{t k}-\delta_t^k$,
$\delta_t^k=\hat \delta_t^k, {{{\underline \Psi}}}_m^k={\hat{{\underline \Psi}}}_m^k,
{{{\overline \Psi}}}_m^k =
{\hat{{\overline \Psi}}}_m^k)$ (because this primal-dual pair satisfies \varepsilonqref{kkt1}). Therefore for $k \geq k_0$
the optimal value ${ z}eta_t^k$ of ${ z}eta_t$
in \varepsilonqref{VFtk} is 0 which gives
${\overline V}_t^{F, k-1}(\pi_n^k)=\mathcal{S}um_{m \in C(n)} p_m\left[b_m^\top \pi_m^k + V_{t+1}^{k-1}(\pi_m^k)\right]$, which, plugged into
gives \varepsilonqref{VFtk}, gives
\begin{equation}gin{equation}
\begin{equation}gin{array}{lcl}
V_t^k(\pi_n^k)-V_t(\pi_n^k)& \leq & \displaystyle
\mathcal{S}um_{m \in C(n)} p_m\left[b_m^\top \pi_m^k + V_{t+1}(\pi_m^k)\right]-V_t(\pi_n^k)\\
&& + \displaystyle  \mathcal{S}um_{m \in C(n)} p_m\left[V_{t+1}^{k-1}(\pi_m^k)-V_{t+1}(\pi_m^k)\right] \\
& \leq & \displaystyle \mathcal{S}um_{m \in C(n)} p_m\left[V_{t+1}^{k-1}(\pi_m^k)-V_{t+1}(\pi_m^k)\right],
\varepsilonnd{array}
\varepsilonnd{equation}
where the last inequality comes from the fact that
$\mathcal{S}um_{m \in C(n)} p_m\left[b_m^\top \pi_m^k + V_{t+1}(\pi_m^k)\right] \leq V_t(\pi_n^k)$
by definition of $V_t$ and feasibility of $\pi_m^k$.

By the induction hypothesis,
$\lim_{k\rightarrow +\infty} V_{t+1}^{k}(\pi_m^k)-V_{t+1}(\pi_m^k)=0$ for all $m \in C(n)$
which, using the equicontinuity of $(V_{t+1}^k)_k$ and the
compacity of $[\underline{\pi}_t, \overline \pi_t]$
implies $\lim_{k\rightarrow +\infty} V_{t+1}^{k-1}(\pi_m^k)-V_{t+1}(\pi_m^k)=0$ for all $m \in C(n)$.
Together with \varepsilonqref{proofconvdsddp}
this implies
$\lim_{k\rightarrow +\infty, k \in \mathcal{S}_n} V_t^k(\pi_n^k)-V_t(\pi_n^k)=0$
and achieves the proof.
$
 \mathcal{S}quare$}\fi

\varepsilonnd{document}