Yan and Reiffers-Masson

CEC-Based Heuristics in Stochastic Optimization Problems

Certainty Equivalence Control-Based Heuristics in Multi-Stage Convex Stochastic Optimization Problems

Chen Yan [Corresponding author]
STATIFY, Inria, 38334 Saint Ismier, France; Biostatistics and Spatial Processes, INRAE, 84914 Avignon, France<EMAIL_ADDRESS>Alexandre Reiffers-Masson
IMT Atlantique, Lab-STICC, UMR CNRS 6285, Brest, France<EMAIL_ADDRESS>
We examine a multi-stage stochastic optimization problem characterized by stagewise-independent, decision-dependent noises with strict constraints. The problem assumes convexity in that, following a specific relaxation, it transforms into a deterministic convex program. The relaxation process is inspired by the principle of Certainty Equivalence Control, which substitutes uncertainties with their nominal values and requires the hard constraints to be satisfied only in an expected sense. Utilizing the solutions obtained from these convex programs, we propose two universal methodologies—re-solving-based and projection-based—to formulate feasible policies relevant to the original problem. These methodologies are subsequently amalgamated to develop a hybrid policy, equipped with a tuning parameter that manages the frequency of re-solving. We derive upper bounds on the gap between the performance of these heuristic policies and the optimal one. Under the Lipschitz-type regularity of the model, these bounds are proportional to the square root of the stochastic noise variance. Assuming additional $\calC^2$-smoothness regularity, an alternative bound, proportional to the variance of the stochastic noise, can be established—providing a refinement when variances are small. Our model provides a framework for dynamic decision-making under uncertainty, encompassing classic inventory and Markovian bandit problems while embracing a broader range of stochastic optimization challenges. We demonstrate our methods using numerical experiments on a network utility maximization problem.

multi-stage stochastic optimization; convex programming; certainty equivalence control

§ INTRODUCTION

Multi-stage stochastic optimization, with its intricate mix of uncertainty and decision-making across time, serves as a cornerstone in a multitude of real-world applications, including supply chain management, power systems operation, inventory control, and financial planning ([Shapiro et al., 2021, Pflug  Pichler, 2014, Küchler, 2009]). This paper focuses on a distinct variant of this type of problem, distinguished by convexity, stagewise-independent decision-dependent noises, and stringent constraints.

To tackle the inherent complexity of the problem, we utilize a specific relaxation process that transforms it into a deterministic convex program. This process is twofold. The first part, inspired by the well-established Certainty Equivalence Control (CEC) principle, replaces uncertainties with their nominal values—typically the mean—and acts as if these were the actual values. It has been demonstrated that, under certain conditions for the linear-quadratic stochastic control problem, this approach actually results in optimal control ([Simon, 1956, Theil, 1957, Duchan, 1974, Runggaldier, 1981]). The second part entails relaxing stringent constraints to be met only in an expected sense. This removes all uncertainties from the equation and offers a valuable approximation in handling intricate stochastic dynamics.

The potency of the approach discussed above, within a multi-stage optimization framework, has been well demonstrated in inventory management ([Kunnumkal  Topaloglu, 2008, Jasin  Kumar, 2012, Cooper, 2002, Secomandi, 2008, Bumpensanti  Wang, 2020]) and Markovian bandits ([Hu  Frazier, 2017, Zayas-Cabán et al., 2017, Brown  Smith, 2020, Zhang  Frazier, 2021, Gast et al., 2021]) contexts, where the relaxed problems emerge as linear programs. Building on these successful applications, we broaden this methodology to encompass a more extensive range of problems, where the relaxed problem can now be a convex program. We propose two universal methodologies—re-solving and projection—to devise feasible policies applicable to the original problem. We then integrate these methodologies to derive a hybrid policy, featuring a tuning parameter to effectively regulate the frequency of re-solving.

In quantifying the effectiveness of these heuristic policies, we establish upper bounds on the gap between their performance and the optimal one. Establishing these upper bounds involves a meticulous manipulation of the propagation of stochastic errors, supplemented with concentration-type inequalities. Assuming Lipschitz-type regularity in the model, these bounds are proportional to the square root of the stochastic noise variance—a first-order estimation. An alternative bound, proportional to the variance of the stochastic noise, is established under the additional $\calC^2$-smoothness regularity. This second-order estimation provides a refinement when variances are small. Such regularity conditions are met on the model, given that the relaxed programs satisfy specific constraint qualifications. This analysis is facilitated by tools from parameterized optimization ([Bonnans  Shapiro, 2013] and <cit.>).

§.§ Related Works and Our Contributions

Our problem assumes an inherent convexity, demonstrated by the CEC-based relaxation process, transforming it into a deterministic convex program. As per our knowledge, previous applications of CEC in multi-stage stochastic optimization settings, such as inventory management or Markovian bandits mentioned previously, restrict themselves to linear programs. The inherent linearity allows for scaling within the model, with the scaling parameter governed by the system size and facilitating asymptotic limit analysis when the size grows large. In a more generalized convex setting, however, such scaling ceases to be applicable in general. Therefore, we shift the perspective to consider a single stochastic problem, presuming that the amplitude of the variance of stochasticities can be reduced—an approach similar to a density model in the linear case (see Remark <ref> for details). Unquestionably, the absence of linearity also leads to more complex analysis and calculations.

We posit that this convexity requirement, to some extent, represents the most general framework to consider under this methodology. Several factors underpin this claim. Primarily, Jensen's inequality applies, ensuring the value of the relaxed convex program is larger than the value of the optimal policy, which is typically elusive. Without this relationship, the upper bounds on the performance gap developed in this paper would become untenable (see Remark <ref> for details). The second factor is computational: convex programs straddle the boundary between efficiently solvable and intractable problems <cit.>. The appeal of applying CEC diminishes if significant challenges already arise at the level of solving the relaxed problems.

In our model, the constraints are "hard": we cannot tolerate any violations, even minor ones. A major challenge we face is that the relaxed problem merely approximates the actual problem, and the decisions derived from it are generally not even feasible for the original problem. Thus, devising simple, efficient, and high-performing feasible policies for the original problem based on the solutions to the relaxed problems is a central theme in these CEC-based techniques. Examples include the LP-index policy ([Gast et al., 2021]) and randomized activation control policy ([Zayas-Cabán et al., 2017]) for the Markovian bandit problem, as well as the booking limit policy and bid-price policy for inventory management ([Jasin  Kumar, 2013]). Another basic idea, based on re-solving and originating from Model Predictive Control ([Rakovic  Levine, 2018]), is also well-documented ([Bumpensanti  Wang, 2020, Gast et al., 2022, Brown  Zhang, 2022]). In this paper, we extend the re-solving concept, termed the "update policy," to a more general convex setting. Simultaneously, we also introduce, to the best of our knowledge, a new concept of projecting the relaxed solution onto the feasible set, which we term the "projection policy", thereby contributing another method for such constructions.

A noteworthy feature of our model, which sets it apart from classical multi-stage stochastic optimization problems as in [Shapiro et al., 2021, Pflug  Pichler, 2014, Küchler, 2009], is the introduction of decision-dependent noise. This divergence is not purely academic; instead, it is motivated by a range of theoretical and practical applications we aim to address—from Markov decision processes ([Puterman, 1994]), stochastic approximation ([Kushner  Yin, 2003]), computational complexity ([Papadimitriou, 1985]), to reinforcement learning ([Sutton  Barto, 2018]). The decision-dependent noise renders various techniques used to handle a standard multi-stage problem, such as quasi-Monte-Carlo approximations and scenario trees, not directly applicable. However, an advantage of the CEC-based heuristics lies in their ability to treat decision-dependent noise in the same manner as i.i.d. white noise, under the assumption of uniformly bounded variances, which we adopt in this paper (see (<ref>)).

Our second-order estimations hinge on local $\calC^2$-smoothness regularity. This characteristic has previously been noted in the simpler linear program case, as exemplified in [Jasin  Kumar, 2012, Wu et al., 2015, Bumpensanti  Wang, 2020] for the inventory model, and in [Zhang  Frazier, 2021, Gast et al., 2021, Brown  Zhang, 2022, Gast et al., 2023] for the Markovian bandit model. These studies introduce this additional regularity under various names such as non-degeneracy or non-singularity, with their precise definitions varying across different papers and contexts. We argue that this regularity property can be expressed as generally as in any finite-dimensional variational inequality ([Facchinei  Pang, 2003]), ultimately establishing a diffeomorphism that renders the constrained parameterized optimization problem locally akin to an unconstrained, parameter-free problem within a neighborhood of the origin in an appropriately dimensioned Euclidean space. We provide a unifying discussion in Appendix appendix:policy-mappingEC.3, building on existing works such as [Robinson, 1982, Robinson, 1987, Dunn, 1987]. We anticipate this unified viewpoint may lend further insight into this property, especially considering its ability to accelerate convergence rates in the asymptotic regime (see Section <ref> for details).

*Outline  The rest of the paper unfolds as follows. We introduce the general model in Section <ref>, and subsequently present three motivating examples in Section <ref>: network utility maximization, inventory management, and Markovian bandit. The concept of CEC is articulated for the general model in Section <ref>. Section <ref> introduces heuristic policies based on CEC, encompassing the update policy, the projection policy, and a hybrid policy featuring a tuning parameter. We present upper bounds on the performance gap of the update and projection policies, predicated on specific regularity conditions, in Section <ref>. Numerical experiments focused on the network utility maximization example are offered in Section <ref>. We conclude with additional comments and a discussion of future research directions in Section <ref>. Proofs, extensions, and additional discussions are in the electronic companion appendix.

*Notational Convention

We use bold letters to denote vectors and matrices, and vectors are represented as row vectors. Capital letters are used to denote random quantities and lowercase letters are used to denote deterministic quantities. The letters $\bx,\bu, \bW, \calE$ are reserved to represent system state, system control, exogenous and endogenous uncertainties, respectively. A norm $\norme{\cdot}$ without subscript is understood to be the $\calL^2$-norm for a vector, and the spectral norm for a matrix, so that it is compatible with the $\calL^2$-norm for vectors in the context of matrix-vector multiplication. We use $(\cdot)^\top$ to denote the transpose, so for instance $\bzero$ is a zero row vector, and $\bonet$ is a column vector of one's. We adopt the convention that time-step begins at $t=1$. For $1 \le t \le T$, denote by $\calF_t$ the $\sigma$-algebra  generated by the random quantities up to time-step $t$, prior to the system state transition, with the convention that $\calF_0 = \{ \bx \}$, where $\bx$ is the initial configuration of the stochastic system (the filtration $\calF_t$ is given more precisely in Assumption <ref> below). We write $\mathbb{E}_{\mathrm{update}} [\cdot]$ (resp. $\mathbb{E}_{\mathrm{proj}} [\cdot]$) to mean the expectation taken under the update policy (resp. the projection policy) that we shall analyse in this work. Suppose that $\bx$ takes values in a domain $\calX$, we denote by $\calB(\bx,\varepsilon)\defeq\left\{\bx' \in \calX \mid \norme{\bx' - \bx} \le \varepsilon \right\}$ the neighbourhood of $\bx$ of radius $\varepsilon$ in $\calX$. Denote by $(\bx,\bw,\bu)$ the row vector that concatenates three row vectors $\bx$, $\bw$ and $\bu$. We write $\bx[t,T]$ for the sequence of vectors $\bx(t')$ for $t \le t' \le T$, so that $\bx[t,T]$ is a short hand for the concatenation of vectors $\big( \bx(t), \bx(t+1), \dots, \bx(T)\big)$. For constraints that involve random quantities, we write a.s. to mean that they are to be satisfied almost surely.

§ THE GENERAL MODEL

Consider the following $T$-stage convex stochastic optimization problem with hard constraints:

Known parameters: initial condition of the system $\bx(1)$; horizon $T$; random vector $\bW$ with known distribution function $f(\bw)$; constraint functions $g_{t,i}(\cdot), h_{t,j}(\cdot)$ for $1 \le t \le T$ and $1 \le i \le I(t)$, $1 \le j \le J(t)$; system evolution with known Markovian laws $\phi(\cdot) + \calE(\cdot)$, where we have separated into two parts for later purpose: the deterministic part $\phi(\cdot)$ is affine, the decision-dependent noise part $\calE(\cdot)$ is with zero mean; reward functions $R_t(\cdot)$ for $1 \le t \le T$.

For each time-step $t=1, \dots, T$:

* The decision-maker gets full knowledge of the current system state $\bX(t)$;

* (The environment) independently draws $\bW(t) \sim f(\bw)$;

* Once $\bW(t)$ has been observed, the decision-maker chooses a control $\bU(t)$ that satisfies the $I(t) + J(t)$ constrains $g_{t,i}(\bX(t), \bW(t), \bU(t)) \le 0 \ a.s.$ and $h_{t,j}(\bX(t), \bW(t), \bU(t)) = 0 \ a.s.$;

* The decision-maker collects a reward $R_t(\bX(t),\bW(t),\bU(t))$;

* The system evolves to the next state $\bX(t+1)$ such that $\bX(t+1) \sim \phi \left( \bX(t), \bW(t), \bU(t) \right) + \calE(\bX(t), \bW(t), \bU(t))$.

Objective: Maximize the expected total sum of rewards over the $T$ time-steps.

Mathematically, the problem can be formulated as follows:
\beginmaxi!|s|[1,T]𝔼 [ ∑_t=1^T R_t ( (t), (t), (t) )  ]  V_opt ((1),T) =
(1) = (1)  a.s.
g_t,i((t), (t), (t)) ≤0   a.s.
h_t,j((t), (t), (t)) = 0   a.s.
(t+1) = ϕ( (t), (t), (t) ) + ((t), (t), (t))  a.s.     1 ≤t ≤T-1
\endmaxi!
where for all $1 \le t \le T$ the model satisfies the following assumptions:

* $\bx(t), \bX(t) \in \R^{n_x}$,  $\bU(t) \in \R^{n_u}$, $\bW(t) \in \R^{n_w}$ are continuous-valued vectors, and are interpreted respectively as the system state (or configuration), system control (or decision, action), and exogenous uncertainties, with $\bx(1)$ being the deterministic initial system configuration. We assume in addition that they all belong to convex subsets of the corresponding Euclidean spaces.

* $R_t: \R^{n_x} \times \R^{n_w} \times \R^{n_u} \rightarrow \R$ are real-valued concave and $\calC^3$-smooth functions for all $t$, jointly for the three arguments. They are interpreted as reward or utility, and are additive across times.

* $\bW(t)$ are bounded independent and identically distributed  (i.i.d.) random vectors with distribution function $f(\bw)$. We write

\begin{equation*}
\expect{\bW} \defeq \barw \in \R^{n_w} \mbox{ and } \var{\bW} \defeq \expect{(\bW - \barw) \cdot (\bW - \barw)^\top}
\end{equation*}

They are interpreted as exogenous uncertainties prior to the decision-making $\bU(t)$ at time-step $t$.

* $g_{t,i}: \R^{n_x} \times \R^{n_w} \times \R^{n_u} \rightarrow \R$ are real-valued convex and $\calC^2$-smooth functions for all $t$ and $1 \le i \le I(t)$, jointly for the three arguments, and $h_{t,j}: \R^{n_x} \times \R^{n_w} \times \R^{n_u} \rightarrow \R$ are real-valued affine functions for all $t$ and $1 \le j \le J(t)$, jointly for the three arguments. They are interpreted as hard constraints to the system and $I(t), J(t) \in \mathbb{N}$ are the numbers of (in)equality constraints at time-step $t$.

* Given $\bX(t), \bW(t)$ and $\bU(t)$, the evolution from $\bX(t)$ to $\bX(t+1)$ is Markovian with an affine behavior in expectation. Consequently, denote by $\calF_t$ the $\sigma$-algebra generated by the random vectors $\bX(t'), \bW(t')$, $\bU(t')$ for $1 \le t' \le t$, plus the random vectors $\calE(\bX(t'), \bW(t'), \bU(t'))$ for $1 \le t' \le t-1$ (note that $\calE(\bX(t), \bW(t), \bU(t))$ is not included in $\calF_t$, but rather in $\calF_{t+1}$),  then

\begin{equation}\label{eq:law-of-Markovian-evolution}
\expect{\bX(t+1) \mid \calF_t} = \expect{\bX(t+1) \mid \bX(t), \bW(t), \bU(t)} \defeq \phi \left( \bX(t), \bW(t), \bU(t) \right)
\end{equation}

where $\phi: \R^{n_x} \times \R^{n_w} \times \R^{n_u} \rightarrow \R^{n_x}$ is an affine function: there exists $\bC \in \R^{(n_x+n_w+n_u) \times n_x}$ and $\bD \in \R^{n_x}$ such that

\begin{equation}\label{eq:phi-explicit}
\phi \left( \bX(t), \bW(t), \bU(t) \right) = (\bX(t),\bW(t),\bU(t)) \cdot \bC + \bD
\end{equation}

A priori $\bC$ and $\bD$ may also depend on $t$. The random decision vector $\bU(t)$ is $\calF_t$-measurable, and therefore cannot depend on any information beyond that available at time-step $t$, adhering to the principle of non-anticipative constraints.

* Write $\calE(\bX(t), \bW(t), \bU(t)) \defeq \bX(t+1) - \phi \left( \bX(t), \bW(t), \bU(t) \right)$, we deduce from Assumption <ref> that

\begin{align*}
& \expect{\calE(\bX(t), \bW(t), \bU(t)) \mid \calF_t} = \expect{\calE(\bX(t), \bW(t), \bU(t)) \mid \bX(t), \bW(t), \bU(t)} \\
& = \int_{\bx \in \R^{n_x}} \bx \ d \nu( \bx \mid \bX(t), \bW(t), \bU(t)) = \mathbf{0}
\end{align*}

where $\nu( \cdot \mid \bX(t), \bW(t), \bU(t))$ is a probability distribution on $\R^{n_x}$ parameterized by $(\bX(t), \bW(t), \bU(t))$. Denote furthermore by

\begin{align*}
& \var{\calE(\bX(t), \bW(t), \bU(t))} = \var{\calE(\bX(t), \bW(t), \bU(t)) \mid \bX(t), \bW(t), \bU(t)} \\
& = \int_{\bx \in \R^{n_x}} \bx \cdot \bx^\top \ d \nu(\bx \mid \bX(t), \bW(t), \bU(t))
\end{align*}

We suppose that

* The support of $\nu(\cdot \mid \bx,\bw,\bu)$ is bounded uniformly on $(\bx,\bw,\bu)$

which implies in particular that

\begin{equation} \label{eq:law-of-stochastic-noise}
\var{\calE} \defeq \sup_{\bx, \bw, \bu} \var{\calE(\bx, \bw, \bu)} < \infty
\end{equation}

$\calE(\bX(t), \bW(t), \bU(t))$ are interpreted as endogenous uncertainties posterior to the decision-making $\bU(t)$ at time-step $t$. They are decision-dependent noises.

* The problem has relatively complete recourse <cit.>. More precisely, for any system configuration $\bX(t)$, any realisation $\bW(t)$, there exists at least one decision $\bU(t)$ that satisfies the hard constraints (<ref>) and (<ref>) at time-step $t$. This assumption is to ensure that the constrained problem (<ref>) is feasible at any moment, so that the feasible region $\calU_t(\bx, \bw)$ given in (<ref>) are always non-empty.

In the above formulation, for each time-step $t$, just before taking our decision, we observe a realization of the random vector $\bW(t)$. The vectors $\bX(t),\bW(t),\bU(t)$ are constrained by $I(t)$ inequality convex functions given by (<ref>), plus $J(t)$ equality affine functions given by (<ref>). The system then evolves to the next state $\bX(t+1)$ according to (<ref>) in a Markovian way, which is an affine function of $(\bX(t),\bW(t),\bU(t))$ for the deterministic part, plus the zero-mean stochastic part $\calE(\bX(t), \bW(t), \bU(t))$ governed by the probability measure $\nu(\cdot \mid \bX(t), \bW(t), \bU(t))$. Our goal is to maximize the sum of rewards over the whole horizon, where the instantaneous reward $R_t(\bX(t),\bW(t),\bU(t))$ being a utility function is concave.

Denote by

\begin{equation}\label{eq:feasible-region-time-t}
\calU_t(\bx,\bw) \defeq \left\{ \bu \in \R^{n_u} \ \big| \ g_{t,i} (\bx,\bw,\bu) \le 0, \ \ \mbox{ for }  1 \le i \le I(t); \ \ h_{t,j} (\bx,\bw,\bu) = 0, \ \ \mbox{ for } 1 \le j \le J(t) \right\}
\end{equation}

which is the set of feasible controls at time-step $t$, given $\bx = \bX(t)$ and $\bw = \bW(t)$. By our previous assumptions on the model, $\calU_t(\bx,\bw)$ is a non-empty (Assumption <ref>), convex (Assumption <ref>), and compact (Assumptions <ref> and <ref>) set in $\R^{n_u}$. We point out that it is possible to relax the boundedness condition on the distributions of $\bW$ and $\calE$ by requiring that they have a tail distribution that converges to zero exponentially fast at infinity, as e.g. a Gaussian distribution, by combining with a concentration inequality as in Lemma <ref> on all our subsequent analysis. We impose the stronger uniformly bounded condition to simplify the matter, since then we can suppose that all the vectors $(\bX(t),\bW(t),\bU(t))$ for $1 \le t \le T$ take values on a common bounded and closed set, and invoke results using compactness.

In order to justify that the optimization problem (<ref>) is well formulated under such generality, we need to impose two additional technical assumptions on the model.

* (Slater CQ for feasible sets) For all $1 \le t \le T$ and all possible realization $(\bx,\bw)$ of $(\bX(t),\bW(t))$, there exists $\bu' \in \calU_t(\bx,\bw)$ such that $g_{t,i}(\bx,\bw,\bu') < 0$ for all $1 \le i \le I(t)$.

* (Lipschitz-continuity on $\nu(\cdot \mid \bx,\bw,\bu)$) There exists $c_{\nu} > 0$ such that for all $(\bx_1,\bw_1,\bu_1)$ and $(\bx_2,\bw_2,\bu_2)$, we have

\begin{equation*}
K(\nu(\cdot \mid \bx_1,\bw_1,\bu_1), \nu(\cdot \mid \bx_2,\bw_2,\bu_2)) \le c_{\nu} \cdot \norme{(\bx_1,\bw_1,\bu_1) - (\bx_2,\bw_2,\bu_2)}
\end{equation*}

where $K(\nu_1,\nu_2)$ is the Kantorovich distance between two probability distributions, or equivalently, the Wasserstain $1$-distance via the Kantorovich duality theorem <cit.>.

We remark that the Slater CQ is a common assumption in convex programming, and the Lipschitz property on the probability measures is a control on the speed of change of the Markovian transition laws as a function of the input $(\bx,\bw,\bu)$. Under this additional Assumption <ref>, we can justify that the optimization problem (<ref>) is well-formulated, the proof is provided in Appendix appendix:proof-well-definedEC.1.

Under Assumptions <ref>-<ref>, the optimization problem (<ref>) is well-defined and an optimal solution exists. Moreover, the mapping $\bx \mapsto V_{\mathrm{opt}} (\bx,T)$ is a continuous function of $\bx$.

Note that an optimal solution to (<ref>), denoted as $U_{\mathrm{opt}}[1,T]$, is a series of random variables, where $U_{\mathrm{opt}}(t) \in \calF_t$, for $1 \le t \le T$.

The convexity assumption, which plays an essential role in our method, will be justified later in Remark <ref>, after introducing the certainty equivalence control on the problem. Formulated as a multistage stochastic optimization problem, the model has stagewise-independent noise <cit.>. It incorporates two sources of stochasticity, $\bW$ and $\calE$, that have distinct nature and are motivated from different application scenarios, as we shall see in the examples displayed in Section <ref>. Clearly the i.i.d. noise $\bW$ is a special and simplified case from the decision-dependent noise $\calE$. Throughout the analysis of the paper we are only concerned with their first and second moments, and it turns out that this subtlety does not play an essential role in our method. Additinally, if the model lacks decision-dependent noise \(\calE\), the Sample Average Approximation, a Monte Carlo method, can be utilized to approximate the optimal solution of (<ref>) <cit.>. Let us also note that state augmentation can be applied by defining \(\widetilde{\bX} \defeq (\bX,\bW)\). Under this setup, the \(\bW\)-part acts as an exogenous and uncontrollable segment of the system state \(\widetilde{\bX}\), as seen in, for example, [Brown  Zhang, 2022]. While this alternative modeling does not introduce technical changes to our subsequent discussions, separating \(\bX\) and \(\bW\) allows algorithms to operate over a reduced state space <cit.>.

In terms of computational complexity, already for linear two-stage stochastic programs with fixed recourse, it is shown in [Hanasusanto et al., 2016] that it is #P-hard to find an approximate solution with sufficiently high accuracy. In the same paper, it is argued that problems with non-relatively complete recourse is even more challenging to solve <cit.>. Intuitively, this is caused by additional implicit constraints at each time-step to ensure feasibility in the future. One common approach to dealt with this issue is to incorporate a high cost for any violation of the hard constraints. Quite often this can transform the program into an equivalent problem with relatively complete recourse, but instead the reward function will have a huge Lipschitz constant <cit.>. To avoid technical difficulties we hence opt to assume that our problem has relatively complete recourse.

For a general $T$-stage program with stagewise-independent and decision-dependent noise, it is shown in [Dyer  Stougie, 2006] that the problem is PSPACE-hard, by treating $T$ as an input parameter. The statue for the complexity of the simpler problem with decision-independent noise is still open (e.g. with only the $\bW$-part in our model (<ref>)). It is conjectured in [Dyer  Stougie, 2006] that it is PSPACE-hard as well.

§ CONCRETE EXAMPLES

The problem (<ref>) covers several classes of models widely studied in the literature as particular cases. We list three below: The first requires the full modeling generality of (<ref>) and will be studied numerically later in Section <ref>; the two later are with affine rewards and polyhedron constraint sets, originate respectively from inventory management and weakly coupled Markov decision processes.

§.§ Network Resource Allocation and Utility Maximization

In a network utility maximization model ([Shakkottai et al., 2008, Srikant  Başar, 2004, Palomar  Chiang, 2006]), $\bX(t)$ represents the bandwidth occupation of $n_x$ routing paths in a communication network with a known topology, $\bW(t)$ represents the arrival of service demands among the $n_x$ paths, and $\bU(t)$ represents the amount of bandwidth that we allocate to each path. Each of the $n_x$ paths use a certain collection of links in the network, and (<ref>)-(<ref>) refers to the constraints on the link capacities, as well as the delays suffered on each routing path, which are typically non-linear convex functions. (<ref>) describes the stochastic evolution of the bandwidth dynamics. Our goal is to maximize a reward function of transferring data flows via the network. One typical choice is $R_t(\bX) = \sum_{i=1}^{n_x} X_i^{1-\alpha}/(1-\alpha)$ for $\alpha > 0$, called the $\alpha$-fairness utility ([Mo  Walrand, 2000]). This is a dynamic extension of the classical network utility maximization problem.

An example of diamond-shaped network with $3$ paths and $5$ links.

Typically, consider a diamond-shaped routing network that can be represented as a directed graph in Figure <ref>. Each directed edge of the graph is called a link, enumerated by $1 \le l \le 5$. Each link has a maximum bandwidth capacity, denoted as $c_l > 0$. At each discrete time-step $1 \le t \le T$, a certain quantity of demand of bandwidth arrives at the source, which will occupy a specific path to reach its destination. A path is a sequence of consecutive directed links that connect the source to the destination, enumerated by $1 \le p \le 3$. We assume that the allocation of bandwidth is instantaneous, meaning that it immediately occupies all the subsequent links of this path. Denote by $(X_1(t),X_2(t),X_3(t))$ the bandwidth occupation of the three routing paths just arriving at time-step $t$, by $(W_1(t),W_2(t),W_3(t))$ the arrivals of new demands of bandwidth on each path at time-step $t$, and the decision is the allocation $(U_1(t),U_2(t),U_3(t))$ of bandwidth among each path. We suppose that the demand is elastic, so that the non-satisfied demand incurs no cost. A first set of constraints on the model can then be expressed as

\begin{equation}\label{eq:network-model-0}
0 \le U_1(t) \le W_1(t), \ \ \  0 \le U_2(t) \le W_2(t), \ \ \ 0 \le U_3(t) \le W_3(t) \ \ \mbox{ for $1 \le t \le T$}
\end{equation}

\begin{align}
& Y_1(t)\defeq U_1(t) + X_1(t) + U_2(t) + X_2(t) \le c_1, \ \ \  Y_2(t) \defeq U_3(t) + X_3(t) \le c_2  \ \ Y_3(t) \defeq U_2(t) + X_2(t) \le c_3 \nonumber  \\
& Y_4(t) \defeq U_1(t) + X_1(t) \le c_4, \ \ \ Y_5(t) \defeq U_2(t) + X_2(t) + U_3(t) + X_3(t) \le c_5  \mbox{ for $1 \le t \le T$}  \label{eq:network-model-1}
\end{align}

To model the stochasticity, we introduce the notation $\calTN(\mu,\sigma^2,a,b)$ to denote a truncated normal distribution, which represents a normal distribution $\calN(\mu,\sigma^2)$ truncated at the interval $[a,b]$. We assume that for $1 \le p \le 3$ and $1 \le t \le T$, the occupation of bandwidth of path $p$ follows the dynamic

\begin{align} \label{eq:network-model-2}
& X_p(t+1) = (X_p(t) + U_p(t)) \cdot q_p  \nonumber \\
& + \calTN \left( 0, \sigma^2, -(X_p(t) + U_p(t))\cdot \min \{q_p, 1-q_p \}, (X_p(t) + U_p(t)) \cdot \min \{q_p, 1-q_p \} \right)
\end{align}

with $0 < q_p < 1$ and $\sigma > 0$. This modeling approach seeks to replicate the scenario in which the bandwidth of a path $p$ is consumed by numerous discrete flows. Each of these flows has a lifespan, which we model using a geometric distribution with parameter $q_p$. After its lifespan ends, the flow releases its occupation of the links corresponding to its path. These flows operate independently of one another. Equation (<ref>) can be interpreted as a continuous approximation of this dynamic. For simplicity, we have chosen to truncate on a symmetric interval around $0$, ensuring that the truncated normal distribution has a mean of zero. In a similar vein, we suppose that the arrivals of new demands follow the distribution

\begin{equation} \label{eq:network-model-3}
W_p(t) = \bar{w}_p + \calTN \left(0,\sigma^2, -\bar{w}_p, \bar{w}_p \right) \mbox{ for $1 \le p \le 3$ and $1 \le t \le T$}
\end{equation}

with $\bar{w}_p > 0$. Note that all variances of the stochastic uncertainties in the model are governed by a single parameter $\sigma$, which provides convenience for our later numerical study.

Following [Parag et al., 2011], we next introduce the link quality degradation function $d_l(y)$ for a link $l$ with a occupation bandwidth $0 \le y < c_l$ as

\begin{equation} \label{eq:network-model-4}
d_l(y) \defeq \frac{y}{c_l(c_l - y)} = \frac{1}{c_l - y} - \frac{1}{c_l} \ \ \mbox{ for $1 \le l \le 5$ and $0 \le y < c_l$}
\end{equation}

Note that $d_l(\cdot)$ is non-negative, convex, and increases from $0$ to $\infty$ when $y$ ranges in $[0, c_l)$. Other choices of $d_l(y)$, as well as its practical implication are discussed in detail in [Parag et al., 2011]. Note that (<ref>) is motivated by seeing link $l$ as a $M/M/1$ queue with arrival rate $y$ and service rate $c_l$, since then $d_l(y)$ is the expected waiting time in the queue.

We assume that the total degradation experienced along path $p$ is the sum of degradation of links that it traverses, and there is a pre-specified positive value $D_p$ being the maximum degradation that can be tolerated on path $p$ at each time-step. This gives rise to a second set of constraints on the model as

\begin{align}  \label{eq:network-model-5}
& d_1(Y_1(t)) + d_4(Y_4(t)) \le D_1, \ \ \ d_1(Y_1(t)) + d_3(Y_3(t)) + d_5(Y_5(t)) \le D_2, \ \ \ d_2(Y_2(t)) + D_5(Y_5(t)) \le D_3  \nonumber \\
& \mbox{ for $1 \le t \le T$}
\end{align}

The decision-maker aims to maximize the following $\alpha$-fairness utility (with $\alpha > 0$)

\begin{equation}  \label{eq:network-model-6}
\expect{\sum_{t=1}^{T} \sum_{p=1}^{3} \frac{(X_p(t)+U_p(t))^{1 - \alpha}}{1 - \alpha} }
\end{equation}

gained by allocating the bandwidth demands over a finite horizon $T$, while respecting the dynamics and constraints described in (<ref>)-(<ref>).

We suppose that initially the network is in a state $\bX(1) \defeq \bx$ such that if we take $\bU(1) \equiv \bzero$, all constraints in (<ref>), (<ref>) and (<ref>) are satisfied with strict inequality. We then remark that by always taking $U_p(t) \equiv \epsilon$ with a fixed $\epsilon > 0$ small enough for all subsequent time-steps $t$, all inequality constraints in (<ref>)-(<ref>) can be satisfied with strict inequality, so Assumption <ref> and Slater CQ of Assumption <ref> are met for this model. Also by our choice of the truncated normal distributions, other conditions in Assumptions <ref>-<ref> concerning the stochastic part of the model are also satisfied.

We shall study numerically one such example using the heuristics developed in this paper in Section <ref>.

§.§ Network Revenue and Inventory Management

In an inventory management scenario, at each time-step $t$, $\bX(t)$ represents the amount of $n_x$ types of resource available in storage, $\bW(t)$ represents the quantity of arrivals of $n_w$ types of (real-valued) customer demands, and $\bU(t)$ consists of two parts: the first part $\bU_1(t)$ represents the amount of additional resource $j$, for $1 \le j \le n_x$, ordered to replenish the storage, at a cost of $c_j$; the second part $\bU_2(t)$ represents the amount of each type $i$, for $1 \le i \le n_w$, of customer demands that we actually serve at time-step $t$, hence $\bU_2(t) \le \bW(t)$. Each type $i$ of customer demand, for $1 \le i \le n_w$, requires a certain combination $a_{ij}$ units of resources $j$ to be served, for $1 \le j \le n_x$, and produce the amount $r_i$ of profit to the decision-maker. While the customer demands that are not served leave the system and are considered as lost with no extra cost. In addition, the resources that are left in the storage incur a holding cost $h_j$ per unit of time, for $1 \le j \le n_x$. So to summarize the reward can be written as

\begin{equation*}
R_t \left( \bX(t), \bW(t), \bU(t) \right) = - \bU_1(t) \cdot \mathbf{c}^\top + \bU_2(t) \cdot \mathbf{r}^\top - \left( \bX(t) + \bU_1(t) - \bU_2(t) \cdot \mathbf{A} \right) \cdot \mathbf{h}^\top
\end{equation*}

where $\mathbf{A}$ is the matrix with entries $a_{ij}$. There is clearly also the constraint $\bX(t) + \bU_1(t) - \bU_2(t) \cdot \mathbf{A} \ge \bzero$. Hence the optimization problem (<ref>) for this inventory management problem can be formulated as
\beginmaxi!|s|[1,T]∑_t=1^T - _1(t) ·𝐜^⊤+ _2(t) ·𝐫^⊤- ( (t) + _1(t) - _2(t) ·𝐀 ) ·𝐡^⊤   V_opt ((1),T) =
(1) = (1)   a.s.
_2(t) - (t) ≤ a.s.
- (t)  ≤ a.s.
- _1(t) ≤,  - _2(t) ≤ a.s.
(t+1) = (t) + _1(t) - _2(t) ·𝐀  a.s.    1 ≤t ≤T-1
\endmaxi!
which is a particular case with linear rewards, linear constraints and without the $\calE$ stochastic part. This is the model that has been considered in <cit.> as a $T$-stage stochastic optimization problem.

Many of its variants have been studied extensively in the literature. Notably, the extension of the model with customer-choice has been considered in [Kunnumkal  Topaloglu, 2008, Bront et al., 2009, Jasin  Kumar, 2012], which incorporate randomness from customer demands' resource consumptions, and can be modeled using the $\calE$-part: $\bU_2(t)$ is hence replaced everywhere by $\bU_2(t) + \calE(\bU_2(t))$, where $\calE(\bU_2(t))$ are noises with distributions depend on $\bU_2(t)$. A simpler model without the possibility of replenishing the storage has been considered in [Cooper, 2002, Secomandi, 2008, Bumpensanti  Wang, 2020], using a re-solving policy similar to the update policy considered in the current work. A thorough discussion for inventory models can be found in [Talluri et al., 2004, Zipkin, 2000].

§.§ Markovian Bandits and Weakly Coupled Markov Decision Processes

A bandit consists of an integer number of symmetric arms, where each arm has $n_x$ states, that evolves as a Markov decision process (MDP). The whole bandit itself can be seen as a large MDP that consists of $N$ statistically identical sub-MDPs (arms), with $N \in \mathbb{N}$ being an integer. The vector $\bX(t)$ represents the number of arms of the bandit being in each of the $n_x$ states, that sums to $N$. An arm can undertake one among a certain number $n_u$ of actions, where each action $j$, for $1 \le j \le n_u$, consumes a certain combination $c_{jk}$, for $1 \le k \le I(t)$ of resources, with a total availability of $\mathbf{b}(t) \in \R^{I(t)}$ resources. In the meantime, an arm in state $i$, for $1 \le i \le n_x$, undertaking action $j$, for $1 \le j \le n_u$, earns a reward $r_{ij}$ for the decision-maker. The vector $\bU(t)$ represents the number of arms undertaking each action, hence it will be more convenient to treat it as a matrix of dimension $n_x \times n_u$, and naturally $\bU(t) \cdot \bonet = \bX(t)$. The resource constraints can be written compactly as $ \bone \cdot \bU(t) \cdot \mathbf{c} \le \mathbf{b}(t)$, with $\mathbf{c}$ being the matrix of entries $c_{jk}$; while the utility collected at time-step $t$ is $\bone \cdot \bU(t) \cdot \mathbf{r}^\top \cdot \bonet$, with $\mathbf{r}$ being the matrix of entries $r_{ij}$. Denote by "$\mathrm{Multinomial}(n,\bp)$" the multinomial distribution with $n$ trials and probability vector of success $\bp$.  For $1 \le i \le n_x$ and $1 \le j \le n_u$, the Markov evolution of $\bX(t)$ can be written as

\begin{align*}
\bX(t+1) & = \sum_{i=1}^{n_x} \sum_{j=1}^{n_u} \mathrm{Multinomial}(U_{ij}(t), \bP^{(j)}_i ) \\
& =  \underbrace{\sum_{j=1}^{n_u} \sum_{j=1}^{n_u} U_{ij}(t) \cdot \bP^{(j)}_i}_{\phi(\bU(t))} + \underbrace{\sum_{i=1}^{n_x} \sum_{j=1}^{n_u} \mathrm{Multinomial}(U_{ij}(t), \bP^{(j)}_i ) - \sum_{j=1}^{n_u} \sum_{j=1}^{n_u} U_{ij}(t) \cdot \bP^{(j)}_i}_{\calE(\bU(t))}
\end{align*}

where $\bP^{(j)}$ for $1 \le j \le n_u$, are in total $n_u$ transition probability matrices of size $n_x \times n_x$, and $\calE(\bU(t))$ are decision-dependent noise obtained by substituting from the multinomial distributions their mean values $\phi(\bU(t))$, the latter being an affine function. So in summary the optimization problem for the Markovian bandit model is
\beginmaxi!|s|[1,T]∑_t=1^T ·(t) ·𝐫^⊤·  V_opt ((1),T) =
(1) = (1)  a.s.
·(t) ·𝐜 - 𝐛(t) ≤ a.s.
- (t) ≤ a.s.
- (t) ·+ (t) =  a.s.
(t+1) = ∑_i=1^n_x ∑_j=1^n_u Multinomial(U_ij(t), ^(j)_i )   a.s.
\endmaxi!
which has linear rewards, linear constraints and without the $\bW$ stochastic part. Note that in the current model $\bX(t)$ and $\bU(t)$ are discrete and integer-valued, and does not fit into the continuous-valued setting we supposed in problem (<ref>). We refer to Remark <ref> for a further discussion on this issue.

This finite horizon Markovian bandit model has been widely studied in the literature, see e.g. [Hu  Frazier, 2017, Zayas-Cabán et al., 2017, Brown  Smith, 2020, Zhang  Frazier, 2021, Gast et al., 2021], as well as its generalization to weakly coupled Markov decision processes [Adelman  Mersereau, 2008, Carpentier et al., 2020, Gast et al., 2022, Brown  Zhang, 2022, Brown  Zhang, 2022]. In [Papadimitriou  Tsitsiklis, 1999], a problem of routing and scheduling in closed queueing networks, called "Network of Queues" has been studied. This problem can be reformulated into the form of (<ref>), by considering each class of jobs as a state, and each server as an action, with an additional action of being idle. The constraints correspond to the classes of jobs that a particular server can serve. By treating $N$ as an input parameter, the problem is shown to be EXP-complete, provided that the horizon $T$ is exponential in $N$, see <cit.> and <cit.>. The simpler "Restless Bandits" problem that has a single server is proven to be PSPACE-complete <cit.>.

§ THE CERTAINTY EQUIVALENT CONTROL (CEC)

The heuristic policies to problem (<ref>) that we shall discuss in this work, i.e. the update policy in Section <ref>, the projection policy in Section <ref>, are inspired from the certainty equivalent control (CEC). The CEC is in general a sub-optimal control that applies at each stage the control that would be optimal if some or all of the uncertain quantities were fixed at their expected values <cit.>. There is an additional difficulty in the current situation to apply CEC, however, as we are facing a problem with hard constraints that depend on both the current system state $\bX(t)$ and disturbance $\bW(t)$, the feasibility of an action is hence of major concern. A key point in the development of the subsequent sections is centered around how to design feasible controls after taking the expectation.

Based on CEC, we apply the following relaxation to the original problem (<ref>): define $\expect{\bX(t)} \defeq \bx(t)$ and $\expect{\bU(t)} \defeq \bu(t)$ where the expectation is taken with the whole trajectory. From the convexity requirement made in Assumption <ref>, these are well-defined system states and controls. By Assumptions <ref>, <ref>, <ref> and Jensen's inequality, we have

\begin{equation} \label{eq:jensen-concave}
\expect{\sum_{t=1}^{T} R_t \left( \bX(t), \bW(t), \bU(t) \right)} \le \sum_{t=1}^{T} R_t \left( \bx(t), \barw, \bu(t) \right).
\end{equation}

\begin{equation}\label{eq:jensen-convex}
\expect{g_{t,i}(\bX(t), \bW(t), \bU(t))} \ge g_{t,i}(\bx(t), \barw, \bu(t)) \ \ \mbox{ for $ 1 \le t \le T$ and $ 1 \le i \le I(t)$}
\end{equation}

\begin{equation}\label{eq:jensen-affine}
\expect{h_{t,j}(\bX(t), \bW(t), \bU(t))} = h_{t,j}(\bx(t), \barw, \bu(t)) \ \ \mbox{ for $ 1 \le t \le T$ and $ 1 \le j \le J(t)$}
\end{equation}

By Assumption <ref>, the expectations can be interchanged with the affine function $\phi(\cdot)$ in (<ref>). By Assumption <ref>, the stochastic parts $\calE(\bX(t), \bW(t), \bU(t))$ are with zero means. All these considerations lead to the following relaxed mathematical program with decision variables $\bu[1,T]$:

\beginmaxi!|s|[1,T]∑_t=1^T R_t ( (t), , (t) )    V_rel- (,T) =
(1) =
g_t,i((t), , (t)) ≤0
h_t,j((t), , (t)) = 0
(t+1) = ϕ( (t), , (t) )      1 ≤t ≤T-1
\endmaxi!
The subscript "$_{\mathrm{vel}-}$" in the notation $V_{\mathrm{rel}-}$ is to make contrast with the later notation $V_{\mathrm{rel}+}$, emphasizing the fact that the relaxed program (<ref>) is solved without the knowledge of $\bW(1)$ at time-step $1$. More generally, for each decision epoch $1 \le t \le T$, we write $V_{\mathrm{rel}-} (\bx,T+1-t)$ to refer to the relaxed CEC problem at time-step $t$ with current system state $\bx(t) = \bx$ and time-span $[t,T]$. Remark that (<ref>) is a deterministic convex program. From (<ref>) and (<ref>), we see that in (<ref>) we are maximizing a larger objective function over a less restrictive feasible region. We hence obtain the key inequality $V_{\mathrm{opt}} (\bx,T) \le V_{\mathrm{rel}-} (\bx,T)$.

The convexity assumptions made in model (<ref>) merit further justification. Specifically, the concavity of the reward functions $R_t(\cdot)$, the convexity of the inequality constraint functions $g_{t,i}(\cdot)$, and the affinity of the equality constraint functions $h_{t,j}(\cdot)$ as well as the system evolution function $\phi(\cdot)$, together ensure the relaxed problem (<ref>) following the CEC remains a convex program. From a computational point of view, there exist solution methods that efficiently solve every convex optimization program satisfying very mild computability restrictions; in contrast, no efficient universal solution methods for non-convex programs are known, and there are strong reasons to expect that no such methods exist, see <cit.> for a thorough discussion. So it is reasonable to remain in the convex paradigm, as otherwise even the relaxed problem may not be tractable.

More importantly, these properties are pivotal as they guarantee that in (<ref>) we are maximizing a larger objective function over a superset of the original feasible region, securing the vital inequality $V_{\mathrm{opt}} (\bx,T) \le V_{\mathrm{rel}-} (\bx,T)$. The significance of maintaining $V_{\mathrm{opt}} (\bx,T) \le V_{\mathrm{rel}-} (\bx,T)$ lies in providing an upper-bound estimate for the sub-optimality gap of any heuristic policy under consideration. This relationship is encapsulated in the inequality:

\begin{align}  \label{eq:pivotal-relation}
\mbox{sub-optimality gap} & \defeq V_{\mathrm{opt}} (\bx,T) - \mbox{value of an heuristic policy} \nonumber \\
& \ \le V_{\mathrm{rel}-} (\bx,T) - \mbox{value of an heuristic policy}
\end{align}

As previously noted in Remark <ref>, determining an exact value for $V_{\mathrm{opt}} (\bx,T)$ to evaluate the sub-optimality gap is generally a complex task. However, $V_{\mathrm{rel}-}(\bx,T)$ can be acquired by solving a convex program, and the value of an efficient heuristic policy can be estimated via Monte-Carlo simulation. Thus, for the remainder of this paper, we will use the right-hand side of (<ref>) as an upper-bound estimate of the sub-optimality gap, termed an optimality gap bound.

§ HEURISTIC POLICIES BASED ON CEC

This section contains the algorithmic results of this paper. Having introduced the CEC in Section <ref>, we construct two heuristic policies based on this principle, the first given in Section <ref> is called the update policy, which requires re-solving a new relaxed mathematical program at each time-step. The second given in Section <ref> is called the projection policy, which solves a single program at the start, and simpler Euclidean projection problems at each later time-step. A hybrid policy that combines the advantages of both of these two policies is then introduced in Section <ref>.

§.§ The Update Policy with Re-Solving

Denote by $\bu^*_{\bx,T}$ an optimal solution of (<ref>), which exists by Proposition <ref>. We use the subscript "$_{\bx,T}$" to keep track that the solution is with respect to initial system configuration $\bx$ and horizon length $T$. The observation is that the first control $\bu^*_{\bx,T}(1)$ from this solution is in general not feasible to the original problem (<ref>), due to the fact that the realization of $\bW(1)$ is still unknown at that moment, and the planning is only guaranteed to be feasible were it be that $\bW(1) = \barw$. The key is that we should apply the re-solving after the realization of $\bW(1)$ is known.

More generally, for each decision epoch $1 \le t \le T$, just before taking a decision, the decision-maker observes the realization $\bw$ of $\bW(t)$ as well as the system state $\bx$, and then solves the following relaxed problem:

\beginmaxi!|s|[t,T]R_t ( (t), , (t) ) + ∑_t'= t+1^T R_t' ( (t'), , (t') )    V̂_rel+ (,T+1-t,) =
(t) =
g_t,i((t), , (t)) ≤0
h_t,j((t), , (t)) = 0
(t+1) = ϕ( (t), , (t) )
g_t',i((t'), , (t')) ≤0
h_t',j((t'), , (t')) ≤0
(t'+1) = ϕ( (t'), , (t') )      t+1 ≤t' ≤T-1
\endmaxi!
Note that the only difference between $V_{\mathrm{rel}-} (\bx,T+1-t)$ and $\hat{V}_{\mathrm{rel}+} (\bx,T+1-t,\bw)$ is that the latter has taken the information $\bW(t) = \bw$ into account, and we actually have

\begin{equation*}
V_{\mathrm{rel}-} (\bx,T+1-t) = \hat{V}_{\mathrm{rel}+} (\bx,T+1-t,\barw)
\end{equation*}

Denote by $\bu^*_{\bx,T+1-t,\bw}[t,T]$ an optimal solution of (<ref>), which exists by Proposition <ref>. We use the subscript "$_{\bx,T+1-t,\bw}$" to keep track that the solution is with respect to system configuration $\bx$, horizon $T+1-t$ and upon observation of $\bW(t)=\bw$. The first control $\bu^*_{\bx,T+1-t,\bw}(t)$ from this control sequence is by construction feasible to (<ref>). Note that we write $\bu^*_{\bx,T+1-t,\bw}(t)$ rather than $\bu^*_{\bx,T+1-t,\bw}(1)$, since although it is the first control to problem (<ref>), it is the $t$-th control to the original problem (<ref>). We repeat this process at each decision epoch $1 \le t \le T$, and this gives rise to the update policy described in Algorithm <ref>.

Initial system configuration vector $\bx(1)$ and time horizon $T$.
Set $\bx \defeq$  current system configuration vector
$t = 1,2,\dots,T$
Observe the realization $\bw$ of $\bW(t)$
Solve $\hat{V}_{\mathrm{rel}+} (\bx,T+1-t,\bw)$ and obtain the control sequence $\bu^*_{\bx,T+1-t,\bw}[t,T]$
Use control $\bu^*_{\bx,T+1-t,\bw}(t)$ to advance to the next time-step
Set $\bx \defeq$  current system configuration vector

The Update Policy.

We use "update" as a short hand for the update policy defined in Algorithm <ref>, and denote by $V_{\mathrm{update}} (\bx,T)$ the value of the update policy with initial condition $\bx$ and horizon $T$.

§.§ The Projection Policy without Re-Solving

The update policy in Algorithm <ref> requires re-solving a convex program that spans the horizon $[t,T]$ at each time-step $t$. This re-solving procedure adjusts the decision based on current available information, and may help to avoid the accumulation of stochastic noises. Yet, from a computational point of view it is not efficient. In this section we propose the idea of policy projection, that replaces the task of re-solving a complex convex program by computing a simpler Euclidean projection.

As before, by solving (<ref>) we obtain an optimal solution $\bu^*[1,T]$ with the corresponding $\bx^*[1,T]$. For each time-step $t$, after the realization of $\bX(t) = \bx(t)$ and $\bW(t) = \bw$, the control $\bu^*(t)$ is in general not feasible, as it was planned with respect to $\bW(t) = \barw$ and $\bX(t) = \bx^*(t)$. In other words, recall the following notation for the set of feasible controls at time-step $t$:

\begin{equation}  \label{eq:feasible-set-for-projection}
\calU_t(\bx,\bw) = \left\{ \bu \in \R^{n_u} \ \big| \ g_{t,i} (\bx,\bw,\bu) \le 0, \ \ \mbox{ for }  1 \le i \le I(t); \ \ h_{t,j} (\bx,\bw,\bu) = 0, \ \ \mbox{ for } 1 \le j \le J(t) \right\}
\end{equation}

then $\bu^*(t) \in \calU_t (\bx^*(t), \barw)$ but in general $\bu^*(t) \notin \calU_t (\bx(t), \bw)$. The idea is that since $\calU_t (\bx(t), \bw)$ is the set of all feasible actions at time-step $t$, which is a non-empty closed set by our assumption, we may apply the Euclidean projection of $\bu^*(t)$ onto $\calU_t (\bx(t), \bw)$ to obtain a feasible action, which is the closest feasible action to $\bu^*(t)$ measured by Euclidean distance. Denote by $\Pi_{\calU_t (\bx(t), \bw)}(\bu^*(t))$ this Euclidean projection. It is the unique solution to the following convex program parameterized by $(\bx(t),\bw)$:

\begin{equation}  \label{eq:convex-program-for-projection}
\min_{\bu} \frac{1}{2} (\bu - \bu^*(t)) \cdot (\bu - \bu^*(t))^\top \mbox{ subject to } \bu \in \calU_t(\bx(t),\bw) \mbox{ given by \eqref{eq:feasible-set-for-projection} }
\end{equation}

We apply this projected action to the system, and repeat this procedure at each time-step $1 \le t \le T$. This gives rise to the projection policy summarized in Algorithm <ref>.

Initial system configuration vector $\bx(1)$ and time horizon $T$.
Solve (<ref>) for time-span $[1,T]$ to obtain an optimal solution $\bu^*[1,T]$
Set $\bx \defeq $  current system configuration vector
$t = 1,2,\dots,T$
Observe the realization $\bw$ of $\bW(t)$
Compute the Euclidean projection $\bu_{\pi}(t) = \Pi_{\calU_t(\bx, \bw)}(\bu^*(t))$
Use control $\bu_{\pi}(t)$ to advance to the next time-step
Set $\bx \defeq$  current system configuration vector

The Projection Policy.

We use "proj" as a short hand for the projection policy defined in Algorithm <ref>, and denote by $V_{\mathrm{proj}} (\bx,T)$ the value of the projection policy with initial condition $\bx$ and horizon $T$.

§.§ A Hybrid Policy

From our previous discussion, we see that the update policy and the projection policy each has its advantage and disadvantage, which are complementary to each other: the update policy is more robust against stochastic uncertainties while requires significantly more computational resource; the projection policy on the other hand, is time-efficient but may not prevent the propagation of estimation errors.

Initial system configuration vector $\bx(1)$ and time horizon $T$. A threshold parameter $\Theta > 0$.
Solve (<ref>) for time-span $[1,T]$ to obtain an optimal solution $\bu^*[1,T]$
Set the current deterministically optimal control as $\mathfrak{U} \defeq \bu^*[1,T]$, and the deterministically optimal system trajectory as $\mathfrak{X} \defeq \bx^*[1,T]$
Set $\bx \defeq$  current system configuration vector
$t = 1,2,\dots,T$
Observe the realization $\bw$ of $\bW(t)$
Compute the Euclidean projection $\bu_{\pi} \defeq \Pi_{\calU_t(\bx, \bw)}(\mathfrak{U}(t))$
Compute the Euclidean distance $\theta \defeq \norme{(\mathfrak{X}(t), \barw, \mathfrak{U}(t)) - (\bx,\bw,\bu_{\pi}) }$
$\theta < \Theta$
Use control $\bu_{\pi}$ to advance to the next time-step
Set $\bx \defeq$ current system configuration vector

Solve $\hat{V}_{\mathrm{rel}+} (\bx,T+1-t,\bw)$ in (<ref>) for time-span $[t,T]$
Update $\mathfrak{U} \defeq \hat{\bu}^*[t,T]$ and $\mathfrak{X} \defeq \hat{\bx}^*[t,T]$ from its solution
Use control $\mathfrak{U}(t)$, which is feasible by construction, to advance to the next time-step
Set $\bx \defeq$ current system configuration vector

A hybrid policy that combines the projection policy and the update policy.

Motivated by this observation, we propose a hybrid algorithm that combines the strength of both policies, which works as follows: we fix a tuning parameter $\Theta > 0$ and keep in the memory a deterministically optimal control $\mathfrak{U}$ and system trajectory $\mathfrak{X}$ from the lastly-solved relaxed convex program. We compute the Euclidean projection as we do in the projection policy. If the "deviation", denoted as $\theta$, of the system to its deterministic counter-part is greater than $\Theta$, then we apply a re-solving as in the update policy. In the meantime we also update $\mathfrak{U}$ and  $\mathfrak{X}$ in the memory. This is summarized in Algorithm <ref>. We remark that as compared to the projection policy, the additional computation of $\theta$ in the hybrid policy is almost free.

Note that in Line 7 of Algorithm <ref>, we have defined $\theta = \norme{(\mathfrak{X}(t), \barw, \mathfrak{U}(t)) - (\bx,\bw,\bu_{\pi})}$ as the deviation at each time-step $t$. We believe that this may capture more information than just using $\norme{\mathfrak{U}(t) - \bu_{\pi}}$. The latter quantity plays an essential role in the theory of error bounds, see e.g. [Luo et al., 1996] and <cit.>, by considering a (non-negative) residual function such as

\begin{equation*}
r_t(\bx, \bw, \bu) \defeq \sum_{i=1}^{I(t)} \left( g_{t,i}(\bx, \bw, \bu)\right)_+ + \sum_{j=1}^{J(t)} \abs{ h_{t,j} (\bx, \bw, \bu) }
\end{equation*}

where $(\cdot)_+ = \max \{ \cdot, 0 \}$. This theory aims at establishing bounds of type

\begin{equation}  \label{eq:resdual-function}
\tau' \cdot r_t(\bx, \bw, \mathfrak{U}(t))^{\gamma'} \le \norme{\mathfrak{U}(t) - \bu_{\pi}} = \norme{\mathfrak{U}(t) - \Pi_{\calU_t(\bx, \bw)}(\mathfrak{U}(t)) } \le \tau \cdot r_t(\bx, \bw, \mathfrak{U}(t))^{\gamma}
\end{equation}

with some positive constants $\tau, \tau'$ and $\gamma, \gamma'$. An immediate application of these error bounds in the implementation of the hybrid policy is the following: the exact value of $\norme{\mathfrak{U}(t) - \bu_{\pi}}$ is not essential, since we only need a comparison of $\theta$ with an empirically chosen threshold $\Theta$. We may very well use the easily computable residual functions in the left and right hand sides of (<ref>) to approximate the quantity $\norme{\mathfrak{U}(t) - \bu_{\pi}}$ for this purpose, provided that these bounds can be proven valid. This saves us from the computation of an Euclidean projection, rendering the hybrid policy even more efficient than the update policy.

§ PERFORMANCE RESULTS

This section contains the major theoretical results of this paper. We state first-order and second-order theorems on the optimality gap bounds for both the update and the projection policies respectively in Sections <ref> and <ref>, and mention briefly the proof ideas. The technical details are given in Appendix appendix:proof-of-main-resultsEC.2. We then discuss interpretation of theses bounds in Section <ref>. Sufficient conditions to satisfy the conditions of these theorems are discussed in length in Appendix appendix:policy-mappingEC.3.

§.§ The Update Policy

Following Remark <ref> and in particular (<ref>), our goal is to compute an upper bound on the difference between $V_{\mathrm{update}} (\bx,T)$ and $V_{\mathrm{rel}-} (\bx,T)$. Our analysis relies heavily on sensitivity analysis of the convex program (<ref>) locally around an optimal solution, which exists by our assumptions. Hence:

\begin{equation}\label{eq:fix-optimal-solution}
\mbox{we fix once and for all an optimal solution $\bu^*[1,T]$ of \eqref{eq:original-problem-rel} with the corresponding $\bx^*[1,T]$}
\end{equation}

For each $1 \le t \le T$, consider the following policy mapping:

\begin{align} \label{eq:policy-mapping}
S^*_t \colon \R^{n_x} \times \R^{n_w} & \rightrightarrows \R^{n_u} \nonumber \\
(\bx,\bw) & \mapsto \left\{ \bu(t) \ \Big| \ \bu[t,T] \mbox{ is an optimal solution to $\hat{V}_{\mathrm{rel}+} (\bx,T+1-t,\bw)$ in \eqref{eq:original-problem-res}}  \right\}
\end{align}

The mapping $S^*_t(\cdot)$ defined above is in general a multi-function, or a set-valued mapping (where "S" stands for solution). We mention that $S_t^*(\cdot)$ are implicit (multi)-functions defined from optimal solutions to a class of parameterized mathematical programs, and in general are not possible to obtain explicit formulas (but see [Gast et al., 2022] and the discussion in Section subsec:non-degeneracy-and-strict-complementarityEC.3.2 for the linear program case). We define the following two well-behaved assumptions on the policy mappings.

[Local Lipschitz-Continuity for $S^*_t(\cdot)$]
For all $1 \le t \le T$, there exists $\varepsilon_t > 0$ such that for all $(\bx,\bw) \in \calB((\bx^*(t),\barw), \varepsilon_t)$, the set $S^*_t(\bx,\bw)$ is single-valued. Moreover, the (single-valued) function $S^*_t(\bx,\bw)$ defined in (<ref>) is locally Lipschitz-continuous in $\calB((\bx^*(t),\barw), \varepsilon_t)$.

[Local $\calC^2$-Smoothness for $S^*_t(\cdot)$]
For all $1 \le t \le T$, there exists $\varepsilon_t > 0$ such that for all $(\bx,\bw) \in \calB((\bx^*(t),\barw), \varepsilon_t)$, the set $S^*_t(\bx,\bw)$ is single-valued. Moreover, the (single-valued) function $S^*_t(\bx,\bw)$ defined in (<ref>) is locally $\calC^2$-smooth in $\calB((\bx^*(t),\barw), \varepsilon_t)$.

We defer the verification and justification of these assumptions to Appendix appendix:policy-mappingEC.3. Specifically, Theorems thm:sufficient-conditions-C2EC.1 and thm:sufficient-conditions-lipschitzEC.2 provide sufficient conditions for satisfying these assumptions. Generally, we anticipate that Assumption <ref> is applicable across a wide variety of contexts. Meanwhile, the more stringent Assumption <ref> depends on additional regularities necessitated by a certain type of Implicit Function Theorem.

We state and prove two performance bound results under each of these assumptions. The key difference between these two results is that in Theorem <ref>, the constant $C_1$ is expressed by the square root of the variances, while in the constant $C_2$ of Theorem <ref>, the square roots have been removed.

Let $V_{\mathrm{opt}} (\bx,T)$ be the value of the stochastic optimization problem (<ref>) that satisfies Assumptions <ref>-<ref>, and let $V_{\mathrm{update}} (\bx,T)$ be the value of the update policy defined in Algorithm <ref>. Under the additional Assumption <ref>, there exists constants $C_1, \mathfrak{P}, \bar{V} > 0$ such that

\begin{equation*}
V_{\mathrm{opt}} (\bx,T) - V_{\mathrm{update}} (\bx,T) \le \mathfrak{P} C_1 +  (1 - \mathfrak{P}) \bar{V}
\end{equation*}

The constant $\mathfrak{P}$ converges to $1$ exponentially fast as both $\var{\bW}$ and $\var{\calE}$ converge to $0$ (an explicit expression is given in (<ref>)). The constant $\bar{V}$ is a finite upper bound of $V_{\mathrm{opt}} (\bx,T)$.
The constant $C_1 = \calO(\sqrt{\var{\bW}} + \sqrt{\var{\calE}})$, with an explicit expression given in (<ref>).

Let $V_{\mathrm{opt}} (\bx,T)$ be the value of the stochastic optimization problem (<ref>) that satisfies Assumptions <ref>-<ref>, and let $V_{\mathrm{update}} (\bx,T)$ be the value of the update policy defined in Algorithm <ref>. Under the additional Assumption <ref>, there exists constants $C_2, \mathfrak{P}, \bar{V} > 0$ such that

\begin{equation*}
V_{\mathrm{opt}} (\bx,T) - V_{\mathrm{update}} (\bx,T) \le \mathfrak{P} C_2 +  (1 - \mathfrak{P}) \bar{V}
\end{equation*}

The constant $\mathfrak{P}$ converges to $1$ exponentially fast as both $\var{\bW}$ and $\var{\calE}$ converge to $0$ (an explicit expression is given in (<ref>)). The constant $\bar{V}$ is a finite upper bound of $V_{\mathrm{opt}} (\bx,T)$.
The constant $C_2 = \calO(\var{\bW} + \var{\calE})$, with an explicit expression given in (<ref>).

Here is the main ingredient of the proofs for these two theorems, details are given in Appendix appendix:proof-of-main-resultsEC.2: we first apply the concentration inequality in Lemma <ref> to bound the probability of the stochastic trajectory remains inside the $\varepsilon_t$-neighbourhood for all $t$, as required in Assumption <ref> (resp. Assumption <ref>). Next, assuming Assumption <ref> (resp. Assumption <ref>), we use the Lipschitz-continuity (resp. $C^2$-smoothness) properties to control the deviations by the variance of the stochastic noises.
The $C^2$-smoothness provides a refinement, since the first-order approximation is linear and the stochastic errors cancel out upon taking expectation, hence only second-order terms remain. See Section <ref> for more discussions.

§.§ The Projection Policy

The analysis of the projection policy relies on understanding how $\Pi_{\calU_t(\bx(t), \bw)}(\bu^*(t))$ behaves for $(\bx(t),\bw)$ in a neighbourhood of $(\bx^*(t),\barw)$, which is a projection onto a perturbed set. Since the Euclidean projection is itself an optimization problem (with a quadratic objective function), we use the same toolkit as in Section <ref>. In abbreviation, we write

\begin{equation}\label{eq:proj}
\proj_t(\bx,\bw) \defeq \Pi_{\calU_t(\bx, \bw)}(\bu^*(t))
\end{equation}

to emphasize the dependence of the projection mapping on $\bx,\bw$ and $t$. Note on the other hand that the vector for applying the projection, $\bu^*(t)$, is fixed for all time-step $t$. Also be definition, $\proj_t(\bx^*(t),\barw) = \bu^*(t)$.

Much like Assumptions <ref> and <ref>, we define the Lipschitz-continuity and the $C^2$-smoothness assumptions on the functions $\proj_t(\cdot)$. A discussion on sufficient conditions to satisfy these assumptions are given in Theorem thm:sufficient-conditions-projectorEC.3, which are specified for the Euclidean projector.

[Lipschitz-Continuity for $\proj_t(\cdot)$]
For all $1 \le t \le T$, there exists $\varepsilon_t > 0$ such that $\proj_t(\bx,\bw)$ as defined in (<ref>) is locally Lipschitz-continuous in $\calB((\bx^*(t),\barw), \varepsilon_t)$.

[$\calC^2$-Smoothness for $\proj_t(\cdot)$]
For all $1 \le t \le T$, there exists $\varepsilon_t > 0$ such that $\proj_t(\bx,\bw)$ as defined in (<ref>) is locally $\calC^2$-smooth in $\calB((\bx^*(t),\barw), \varepsilon_t)$.

The next two theorems are the analogues of Theorem <ref> and Theorem <ref> for the projection policy. Proofs are detailed in Appendix subsec:proof-of-projectionEC.2.2.

Let $V_{\mathrm{opt}} (\bx,T)$ be the value of the stochastic optimization problem (<ref>) that satisfies Assumptions <ref>-<ref>, and let $V_{\mathrm{proj}} (\bx,T)$ be the value of the projection policy defined in Algorithm <ref>. Under the additional Assumption <ref>, there exists constants $C_3, \mathfrak{P}', \bar{V} > 0$ such that

\begin{equation*}
V_{\mathrm{opt}} (\bx,T) - V_{\mathrm{proj}} (\bx,T) \le \mathfrak{P}' C_3 +  (1 - \mathfrak{P}') \bar{V}
\end{equation*}

The constant $\mathfrak{P}'$ converges to $1$ exponentially fast as both $\var{\bW}$ and $\var{\calE}$ converge to $0$ (an explicit expression is given in (<ref>)). The constant $\bar{V}$ is a finite upper bound of $V_{\mathrm{opt}} (\bx,T)$. The constant $C_3 = \calO(\sqrt{\var{\bW}} + \sqrt{\var{\calE}})$, with an explicit expression given in (<ref>).

Let $V_{\mathrm{opt}} (\bx,T)$ be the value of the stochastic optimization problem (<ref>) that satisfies Assumptions <ref>-<ref>, and let $V_{\mathrm{proj}} (\bx,T)$ be the value of the projection policy defined in Algorithm <ref>. Under the additional Assumption <ref>, there exists constants $C_4, \mathfrak{P}', \bar{V} > 0$ such that

\begin{equation*}
V_{\mathrm{opt}} (\bx,T) - V_{\mathrm{proj}} (\bx,T) \le \mathfrak{P}' C_4 +  (1 - \mathfrak{P}') \bar{V}
\end{equation*}

The constant $\mathfrak{P}'$ converges to $1$ exponentially fast as both $\var{\bW}$ and $\var{\calE}$ converge to $0$ (an explicit expression is given in (<ref>)). The constant $\bar{V}$ is a finite upper bound of $V_{\mathrm{opt}} (\bx,T)$. The constant $C_4 = \calO(\var{\bW} + \var{\calE})$, with an explicit expression given in (<ref>).

§.§ Interpretation of the Optimality Gap Bounds

In this subsection, we illustrate via en elementary observation as how the additional $\calC^2$-smoothness can achieve for a refinement, and justify how and when this can be called a refinement via the previously studied examples in Section <ref>.

For this purpose, let $Y$ be a real-valued random variable with mean value $\expect{Y}=y$ and a small variance, with the meaning of "small" that will be made precise in the analysis. Suppose we are interested in finding an upper bound for the quantity $\abs{\expect{v(Y)} - v(y)}$, with a certain real-valued $\calC^2$-smooth function $v(\cdot)$. Let $\bar{v}$ be an upper bound of the function $v(\cdot)$ over the support of $Y$, and we suppose that $\bar{v} < \infty$. Let us choose a small $\varepsilon$-neighbourhood of $y$.  If we only rely on the Lipschitz-continuity of the function $v(\cdot)$, we obtain

\begin{align}
\abs{\expect{v(Y)} - v(y)} & \le L_v \cdot \expect{\abs{Y - y} \cdot \mathds{1}_{Y \in [y-\varepsilon, y+\varepsilon]}} + 2 \bar{v} \cdot \proba{Y \notin [y-\varepsilon, y+\varepsilon]}  \nonumber \\
& \ \le L_v \cdot \varepsilon + C_1 \cdot \exp{ \left( - \frac{C_2 \cdot \varepsilon^2 }{\var{Y}} \right) }  \label{eq:illustration-1}
\end{align}

with $L_v$ being the Lipschitz-constant of $v(\cdot)$ in the interval $[y-\varepsilon, y+\varepsilon]$, and $C_1,C_2 > 0$ are constants from Lemma <ref>. On the other hand, if we use the additional $\calC^2$-smoothness, we have

\begin{equation*}
v(Y) = v(y) + v'(y) \cdot (Y-y) + \frac{1}{2} v^{''}(\tilde{y}_Y) \cdot (Y-y)^2
\end{equation*}

with $\tilde{y}_Y$ being some value depending on $Y$. So we deduce that

\begin{align}
\abs{\expect{v(Y)} - v(y)} & \le \frac{H_v}{2} \cdot \expect{(Y-y)^2 \cdot \mathds{1}_{Y \in [y-\varepsilon, y+\varepsilon]}} + 2 \bar{v} \cdot \proba{Y \notin [y-\varepsilon, y+\varepsilon]} \nonumber \\
& \ \le \frac{H_v}{2} \cdot \varepsilon^2 +  C_1 \cdot \exp{\left( - \frac{C_2 \cdot \varepsilon^2 }{\var{Y}}\right)} \label{eq:illustration-2}
\end{align}

where $H_v$ is some upper found for the function $v''(\cdot)$ in $[y-\varepsilon, y+\varepsilon]$. If $\var{Y}$ is much smaller than $\varepsilon^2$, say $\var{Y} = \varepsilon^3$, then the second term in (<ref>) and (<ref>) is negligible as compared to their respective first terms. If in the meantime, the value of $\varepsilon$ can also chosen to be small enough, the bound in (<ref>) is of first-order $\calO(\varepsilon)$, while the bound in (<ref>) is of second-order $\calO(\varepsilon^2)$. So in summary, a sufficient condition on $\var{Y}$ and $\varepsilon$ to ensure that the upper bound (<ref>) is a refinement to the upper bound (<ref>) is $\var{Y} \ll \varepsilon^2 \ll 1$. Note that this is not a necessary condition and in any case, a second-order bound, if exists, always provides an alternative estimation to the optimality gap. We refer to Figure <ref> for an illustration of this effect in a numerical study of the network utility maximization example.

In Examples <ref> and <ref>, both relaxed problems corresponding to (<ref>) and (<ref>) are linear programs, enabling scaling of the model size, which in turn gives rise to a reduction in variances commensurate with the scaling. For instance, in Example <ref>, the vectors $\bX$ and $\bU$ denote the number of bandit arms in each of the $n_x$ states and taking each of the $n_u$ actions respectively, totaling $N$. Consequently, the vectors $\bMN \defeq \bX/N$ and $\bYN \defeq \bU/N$ represent the corresponding proportions. Given that every (in)equality defining the optimization problem (<ref>) is affine, dividing each formula by $N$, while allowing fractional numbers of arms, transforms it into an equivalent problem involving the proportionality real-valued vectors $\bMN$ and $\bYN$ within some simplex with appropriate dimensions. Specifically, (<ref>) transforms into $\bMN(t+1) = \phi(\bYN(t)) + \bEN(t)$, where $\phi(\cdot)$ represents an affine function and $\bEN(t)$ serves as a random vector constituting the stochastic component, $\calE$. Here, $\bEN(t)$ adheres to $\var{\calE} = \calO(1/N)$, as proven in <cit.>. Thus, as the number of symmetric arms $N$ increases, a second-order result as found in Theorems <ref> and <ref> delivers a faster $\calO(1/N)$ asymptotic convergence rate for the optimality gap bounds, contrasting with the $\calO(1/ \sqrt{N})$ rate seen in Theorem <ref> and Theorem <ref>. In fact, [Gast et al., 2022] has demonstrated that with a perfect rounding scheme to fit the original integer-valued problem (<ref>), the convergence rate can be enhanced to an even faster rate of $\exp(-\calO(N))$.

Similarly, in Example <ref> for the inventory model without the option of replenishing storage, a common scaling method is to multiply both the initial resource and the time horizon by a factor $N$ (not necessarily integer). As shown in [Jasin  Kumar, 2012], a faster $\calO(1/N)$ asymptotic convergence rate can be achieved if the relaxed linear program satisfies a non-degenerate condition. We emphasize that these efficient approximations do not contradict the computational hardness results, e.g. those mentioned in Section <ref>, since the scaling regimes under the two situations are different. For example, the "Restless Bandits" problem in [Papadimitriou  Tsitsiklis, 1999] under this scaling involves identifying the best $5$, $50$, or $500$ arms out of $10$, $100$, or $1000$ arms, respectively. By contrast, identifying a fixed quantity of $5$ arms, regardless of the total number of arms, is PSPACE-complete. Also, note that in the setting of a convex program studied in this paper, such scaling is in general not possible. For situations where the amplitude of the variances are kept constant, thereby making a consideration of robust optimality more suitable, we refer to Remark rem:asymptotic-vs-robustnessEC.1 for further discussion.

§ NUMERICAL RESULTS ON NETWORK UTILITY MAXIMIZATION

Path $p$/Link $l$     x     $\barw$     q     $\alpha$     D     c

1     1     2     0.6     0.5     100     6

2     1     2     0.7     0.5     100     4

3     1     2     0.5     0.5     100     3

4     -     -     -     -     -     4

5     -     -     -     -     -     6

Explicit parameters used in the numerical example of Figure <ref>.

$\Theta \equiv 1.5$

$\Theta = [0.4, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0]$

Fix $\sigma \equiv 1.0$ and vary $\Theta$: the trade-off between policy performance and computation efficiency

The model in Figure <ref> with parameters given in Table <ref>, horizon is $T=30$. The annotated numbers are the counts of re-solving used in the hybrid policy, averaged over $400$ simulations (the initial solving in Line $1$ of Algorithm <ref> before observing $\bW(1)$ is also included). Since $T=30$, this count is upper bounded by $31$.

In this section we provide a numerical study on Example <ref>. The parameters employed in our experiments are detailed in Table <ref>. Our code, written in Python, can be accessed from our GitHub repository (available via <https://gitlab.inria.fr/phdchenyan/network_scheduling.git>). We utilize the CVXPY package for solving convex programs. In addition to the update, projection, and hybrid policies, we also introduce the myopic policy as a benchmark; this policy undertakes a random feasible action at each time-step. This random selection is subject to the specific convex program solver, but is generally considered the least computationally demanding of all feasible policies.

In the first set of experiments, we set the horizon $T=30$ with the aim of comparing the performance of various policies studied in this paper. The simulation results are presented in Figure <ref>. We evaluate the variance variable $\sigma$ within the range of $[0, 2]$. In Figure <ref>, the tuning parameter $\Theta$ for the hybrid policy remains fixed at $\Theta=1.5$, while in Figure <ref>, a different and increasing $\Theta$ is selected for each ascending value of $\sigma$. In Figure <ref>, we hold $\sigma = 1$ constant and adjust $\Theta$ within the scope of [0.5,5]. The annotated numbers indicate the counts of re-solving utilized in each policy, inclusive of the initial solving of the convex program to maintain a uniform comparison. This count is capped at $31$. By definition, it stands at $0$ for the myopic policy, at $1$ for the projection policy, and at $31$ for the update policy. For the hybrid policy, it assumes a value between $1$ and $31$ depending on $\Theta$, with a larger $\Theta$ resulting in fewer re-solving counts.

For the hybrid policy, the re-solving count does not account for when the re-solving occurs, as this is entirely stochastic and reliant on the realization of each sample run. However, it can serve as an indicator of computational resource consumption. From Figure <ref>, we deduce that more frequent re-solving can improve performance, and depending on the variance $\sigma$, a model-driven selection of $\Theta$ can help strike a balance between performance assurance and computational resource consumption. For instance, when $\sigma = 1$, choosing $\Theta = 3$ appears to be a judicious decision for this purpose. Meanwhile, when $\sigma = 0.1$, the projection policy proves ideal: Not only is it substantially more efficient than the update policy, but it also achieves a noteworthy improvement over the myopic policy, merely by solving the initial convex program once. Note that in the asymptotic limit $\sigma = 0$, all policies, except the myopic, achieve performance that coincides with the value of the relaxed convex program.

The model in Figure <ref> with parameters given in Table <ref>, horizon is $T=3$. We plot the optimality gap as function of a single parameter $\sigma$, which quantifies the strength of the variances, from both $\bW$ and $\calE$, for the update policy. We then apply a linear and quadratic fit. The plots are done in two scale, one for $\sigma$ in the range $[0.025, 0.3]$, where a quadratic fit is more suitable; another in the range $[0.2,0.8]$, where a linear fit appears just as adequate as a quadratic fit.

For a second set of experiments, we set the horizon $T=3$ and aim to validate the convergence rate backed by our theoretical results. The results are depicted in Figure <ref>. We begin by checking the sufficient conditions in Theorem <ref> for the update policy. As per Theorem thm:sufficient-conditions-C2EC.1, this boils down to validating LICQ and the strict complementarity to the initial convex program (<ref>). It should be noted that for the current model, the constraints in (<ref>) are actually superfluous and implied by the constraints in (<ref>). Numerically, we discover that for each of the three time-steps, only one constraint in (<ref>) is saturated, with the LICQ being upheld. The strict complementarity is also satisfied, after examining the Lagrange multipliers from the solution to (<ref>). On the other hand, the projection mapping (<ref>) for time-step $t=1$ is degenerate (see Appendix appendix:policy-mappingEC.3 and in particular Figure <ref> for an illustration), suggesting that the rate claimed in Theorem <ref> does not apply to the projection policy.

As a consequence of the preceding analysis, we observe in the left panel of Figure <ref> that for the update policy with $\sigma$ in the range of $[0.025,0.3]$, a quadratic fit of the optimality gap is more suitable than a linear fit. However, in the right panel of Figure <ref>, when we extend the range of $\sigma$ to $[0.2,0.8]$, a linear fit appears just as adequate as a quadratic fit for $\sigma$ in this scope. This aligns with the analysis in Section <ref> and the first step of the proof of Theorem <ref> via concentration-type inequalities: when the variance is substantial, the concentration bound of Lemma <ref> becomes too broad, and the bound in (<ref>) (a first-order linear bound) proves more accurate than the one in (<ref>) (a second-order quadratic bound). We also mention that for the projection policy, the fits appear to be linear for $\sigma$ in both the ranges $[0.025,0.3]$ and $[0.2,0.8]$ (not shown in the figure). This supports the observation that the projection mapping is degenerate at $t=1$, and hence the quadratic convergence rate of Theorem <ref> is not applicable.

§ CONCLUSION AND FUTURE WORKS

We addressed multi-stage stochastic optimization problems that exhibit a convex structure upon applying a CEC-based relaxation. Our framework not only embraces established models from inventory management and Markovian bandits, but also adeptly manages more intricate, non-linear problems like network utility maximization. The heuristic policies we introduced offer a quantifiable performance gap in relation to the optimal ones, while also balance between performance and efficiency. As we conclude, we chart a path for several avenues of future research.

* (Infinite Horizon) In situations where all the data of the model (i.e., rewards, noises, and constraints) are time-independent, we can formulate the problem under the infinite horizon with a time-averaged reward criterion. Specific models of (<ref>), such as the "Restless Bandits", have already been analysed within this framework, as seen with the well-known Whittle index policy ([Whittle, 1988]) and the LP-priority policy ([Verloop, 2016]). A distinguishing feature of the infinite horizon case is that, when appropriately formulated, we only have a single relaxed program that describes the system's stationary behavior - significantly simpler than the finite horizon case, which is size $T$ times larger. In principle, the idea from the projection policy in this paper still applies. However, it now necessitates an additional verification of asymptotic stability to ensure optimality. In a recent paper ([Hong et al., 2023]), which studied the "Restless Bandits" model, the authors circumvented this asymptotic stability condition using a coupling technique. Exploring how to integrate their approach into the more general setting in this paper is a worthwhile endeavor.

* (The Effect of Re-Solving) The numerical experiments in Section <ref> indicate that more frequent re-solving always yields improved performance in this network utility maximization model. Additionally, the analysis in Appendix subsec:Lipschitz-constant-of-V-rel-postEC.4.2 hints at the update policy's potential for increased noise immunity. Nonetheless, a comprehensive understanding of the effects of re-solving is far from trivial, since it requires comparing different approximation schemes for a stochastic problem. Note that this question should not be conflated with the concept of time-consistency in multi-stage stochastic programming <cit.>. Although it may seem counterintuitive, more re-solving is not always beneficial - a fact underscored by counter-examples for the inventory management model found in [Cooper, 2002, Secomandi, 2008, Bumpensanti  Wang, 2020]. We plan to explore this issue further and report our findings in a follow-up paper, utilizing a broader framework than the one presented in this paper (convexity is not needed for studying the effect of re-solving). Our numerical experiments appear to suggest that the negative impact of re-solving is related to degeneracy, as we define in Appendix appendix:policy-mappingEC.3, in conjunction with the presence of stringent constraints.

* (Discrete-Valuedness) One limitation in our model is the assumption of continuous-valued data. However, in many real-world applications, the model is intrinsically discrete-valued, such as when the state is represented by counting numbers as in the Markovian bandit model of Example <ref>, or when the action is binary as in some variants of the inventory management model of Example <ref>. In many cases, continuous-valuedness is more a hindrance than a help, as unless explicit analytical solutions are available, continuous functions must be discretized for numerical solving. In our context, however, it is the inverse: discrete-valuedness disrupts the desired smoothness and convexity properties essential for our CEC-based analysis. Certain situations allow us to overcome these difficulties by relaxing discrete values, e.g. by permitting fractional numbers of arms as discussed in Remark <ref> and [Gast et al., 2021, Gast et al., 2022], or by treating binary action as a probability with a value between $0$ and $1$, as in the probabilistic allocation policy in [Bumpensanti  Wang, 2020].

Chen Yan extends heartfelt gratitude to Nicolas Gast and Bruno Gaujal for their invaluable guidance and numerous enlightening discussions on the intriguing subjects related to this research. He also wishes to express his appreciation to IMT Atlantique (Campus de Brest) for its hospitality. The time spent there facilitated many productive conversations with Alexandre Reiffers-Masson, ultimately making this paper possible.

[Adelman  Mersereau, 2008]
Adelman D, Mersereau AJ (2008) Relaxations of weakly coupled stochastic dynamic
programs. Operations Research 56(3):712–727.

[Avriel  Williams, 1970]
Avriel M, Williams A (1970) The value of information and stochastic
programming. Operations Research 18(5):947–954.

[Ben-Tal et al., 2009]
Ben-Tal A, El Ghaoui L, Nemirovski A (2009) Robust optimization,
volume 28 (Princeton university press).

[Ben-Tal  Nemirovski, 2001]
Ben-Tal A, Nemirovski A (2001) Lectures on modern convex optimization:
analysis, algorithms, and engineering applications (SIAM).

[Bertsekas, 2012]
Bertsekas D (2012) Dynamic programming and optimal control: Volume I,
volume 1 (Athena scientific).

[Bertsimas et al., 2011]
Bertsimas D, Brown DB, Caramanis C (2011) Theory and applications of robust
optimization. SIAM review 53(3):464–501.

[Blondel  Tsitsiklis, 2000]
Blondel VD, Tsitsiklis JN (2000) A survey of computational complexity results
in systems and control. Automatica 36(9):1249–1274.

[Bonnans  Shapiro, 2013]
Bonnans JF, Shapiro A (2013) Perturbation analysis of optimization
problems (Springer Science & Business Media).

[Bront et al., 2009]
Bront JJM, Méndez-Díaz I, Vulcano G (2009) A column generation
algorithm for choice-based network revenue management. Operations
research 57(3):769–784.

[Brown  Smith, 2020]
Brown DB, Smith JE (2020) Index policies and performance bounds for dynamic
selection problems. Manag. Sci. 66:3029–3050.

[Brown  Zhang, 2022]
Brown DB, Zhang J (2022a) Fluid policies, reoptimization, and
performance guarantees in dynamic resource allocation. Reoptimization,
and Performance Guarantees in Dynamic Resource Allocation (November 5, 2022)

[Brown  Zhang, 2022]
Brown DB, Zhang J (2022b) On the strength of relaxations of weakly
coupled stochastic dynamic programs. Operations Research .

[Bumpensanti  Wang, 2020]
Bumpensanti P, Wang H (2020) A re-solving heuristic with uniformly bounded loss
for network revenue management. Management Science 66(7):2993–3009.

[Carpentier et al., 2020]
Carpentier P, Chancelier JP, De Lara M, Pacaud F (2020) Mixed spatial and
temporal decompositions for large-scale multistage stochastic optimization
problems. Journal of Optimization Theory and Applications

[Cooper, 2002]
Cooper WL (2002) Asymptotic behavior of an allocation policy for revenue
management. Operations Research 50(4):720–727.

[Duchan, 1974]
Duchan AI (1974) A clarification and a new proof of the certainty equivalence
theorem. International Economic Review 216–224.

[Dunn, 1987]
Dunn JC (1987) On the convergence of projected gradient processes to singular
critical points. Journal of Optimization Theory and Applications

[Dyer  Stougie, 2006]
Dyer M, Stougie L (2006) Computational complexity of stochastic programming
problems. mathematical programming 106:423–432.

[Facchinei  Pang, 2003]
Facchinei F, Pang JS (2003) Finite-dimensional variational inequalities
and complementarity problems (Springer).

[Gast, 2017]
Gast N (2017) Expected values estimated via mean-field approximation are
1/n-accurate. Proceedings of the ACM on Measurement and Analysis of
Computing Systems 1(1):1–26.

[Gast et al., 2021]
Gast N, Gaujal B, Yan C (2021) Lp-based policies for restless bandits:
necessary and sufficient conditions for (exponentially fast) asymptotic
optimality. arXiv preprint arXiv:2106.10067 .

[Gast et al., 2022]
Gast N, Gaujal B, Yan C (2022) The lp-update policy for weakly coupled markov
decision processes. arXiv preprint arXiv:2211.01961 .

[Gast et al., 2023]
Gast N, Gaujal B, Yan C (2023) Exponential asymptotic optimality of whittle
index policy. Queueing Systems 1–44.

[Hanasusanto et al., 2016]
Hanasusanto GA, Kuhn D, Wiesemann W (2016) A comment on “computational
complexity of stochastic programming problems”. Mathematical
Programming 159:557–569.

[Hong et al., 2023]
Hong Y, Xie Q, Chen Y, Wang W (2023) Restless bandits with average reward:
Breaking the uniform global attractor assumption. arXiv preprint
arXiv:2306.00196 .

[Hu  Frazier, 2017]
Hu W, Frazier P (2017) An asymptotically optimal index policy for
finite-horizon restless bandits. arXiv preprint arXiv:1707.00205 .

[Huang et al., 1977]
Huang CC, Vertinsky I, Ziemba WT (1977) Sharp bounds on the value of perfect
information. Operations Research 25(1):128–139.

[Ioffe, 2017]
Ioffe AD (2017) Variational analysis of regular mappings. Springer
Monographs in Mathematics. Springer, Cham .

[Jasin  Kumar, 2012]
Jasin S, Kumar S (2012) A re-solving heuristic with bounded revenue loss for
network revenue management with customer choice. Mathematics of
Operations Research 37(2):313–345.

[Jasin  Kumar, 2013]
Jasin S, Kumar S (2013) Analysis of deterministic lp-based booking limit and
bid price controls for revenue management. Operations Research

[Kohler  Lucchi, 2017]
Kohler JM, Lucchi A (2017) Sub-sampled cubic regularization for non-convex
optimization. International Conference on Machine Learning, 1895–1904

[Küchler, 2009]
Küchler C (2009) Stability, Approximation, and Decomposition in
Two-and Multistage Stochastic Programming (Springer).

[Kunnumkal  Topaloglu, 2008]
Kunnumkal S, Topaloglu H (2008) A refined deterministic linear program for the
network revenue management problem with customer choice behavior. Naval
Research Logistics (NRL) 55(6):563–580.

[Kushner  Yin, 2003]
Kushner H, Yin G (2003) Stochastic Approximation and Recursive Algorithms
and Applications. Stochastic Modelling and Applied Probability (Springer New
York), ISBN 9780387008943.

[Lu  Robinson, 2008]
Lu S, Robinson SM (2008) Variational inequalities over perturbed polyhedral
convex sets. Mathematics of Operations Research 33(3):689–711.

[Luo et al., 1996]
Luo ZQ, Pang JS, Ralph D (1996) Mathematical programs with equilibrium
constraints (Cambridge University Press).

[Mo  Walrand, 2000]
Mo J, Walrand J (2000) Fair end-to-end window-based congestion control.
IEEE/ACM Transactions on networking 8(5):556–567.

[Palomar  Chiang, 2006]
Palomar DP, Chiang M (2006) A tutorial on decomposition methods for network
utility maximization. IEEE Journal on Selected Areas in Communications

[Papadimitriou, 1985]
Papadimitriou CH (1985) Games against nature. Journal of Computer and
System Sciences 31(2):288–301.

[Papadimitriou  Tsitsiklis, 1999]
Papadimitriou CH, Tsitsiklis JN (1999) The complexity of optimal queuing
network control. Math. Oper. Res 293–305.

[Parag et al., 2011]
Parag P, Sah S, Shakkottai S, Chamberland JF (2011) Value-aware resource
allocation for service guarantees in networks. IEEE Journal on Selected
Areas in Communications 29(5):960–968.

[Pflug  Pichler, 2014]
Pflug GC, Pichler A (2014) Multistage stochastic optimization, volume
1104 (Springer).

[Puterman, 1994]
Puterman ML (1994) Markov Decision Processes: Discrete Stochastic Dynamic
Programming (New York, NY, USA: John Wiley & Sons, Inc.), 1st edition.

[Rakovic  Levine, 2018]
Rakovic SV, Levine WS (2018) Handbook of model predictive control

[Reaiche, 2016]
Reaiche M (2016) A note on sample complexity of multistage stochastic programs.
Operations Research Letters 44(4):430–435.

[Robinson, 1982]
Robinson SM (1982) Generalized equations and their solutions, part II:
applications to nonlinear programming (Springer).

[Robinson, 1987]
Robinson SM (1987) Local structure of feasible sets in nonlinear
programming, Part III: Stability and sensitivity (Springer).

[Robinson, 2003]
Robinson SM (2003) Constraint nondegeneracy in variational analysis.
Mathematics of Operations Research 28(2):201–232.

[Runggaldier, 1981]
Runggaldier WJ (1981) A generalized certainty-equivalence result in stochastic
control. Automatica 17(2):393–395.

[Secomandi, 2008]
Secomandi N (2008) An analysis of the control-algorithm re-solving issue in
inventory and revenue management. Manufacturing & Service Operations
Management 10(3):468–483.

[Shakkottai et al., 2008]
Shakkottai S, Srikant R, et al. (2008) Network optimization and control.
Foundations and Trends® in Networking 2(3):271–379.

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

[Shapiro  Nemirovski, 2005]
Shapiro A, Nemirovski A (2005) On complexity of stochastic programming
problems. Continuous optimization: Current trends and modern
applications 111–146.

[Simon, 1956]
Simon HA (1956) Dynamic programming under uncertainty with a quadratic
criterion function. Econometrica, Journal of the Econometric Society

[Srikant  Başar, 2004]
Srikant R, Başar T (2004) The mathematics of Internet congestion
control (Springer).

[Sutton  Barto, 2018]
Sutton RS, Barto AG (2018) Reinforcement learning: An introduction (MIT

[Talluri et al., 2004]
Talluri KT, Van Ryzin G, Van Ryzin G (2004) The theory and practice of
revenue management, volume 1 (Springer).

[Theil, 1957]
Theil H (1957) A note on certainty equivalence in dynamic planning.
Econometrica: Journal of the Econometric Society 346–349.

[Verloop, 2016]
Verloop M (2016) Asymptotically optimal priority policies for indexable and
nonindexable restless bandits. Annals of Applied Probability

[Villani, 2021]
Villani C (2021) Topics in optimal transportation, volume 58 (American
Mathematical Soc.).

[Whittle, 1988]
Whittle P (1988) Restless bandits: activity allocation in a changing world.
Journal of Applied Probability 25A:287–298.

[Wu et al., 2015]
Wu H, Srikant R, Liu X, Jiang C (2015) Algorithms with logarithmic or sublinear
regret for constrained contextual bandits. Advances in Neural
Information Processing Systems 28.

[Yan, 2022]
Yan C (2022) Close-to-opimal policies for Markovian bandits. Ph.D.
thesis, Université Grenoble Alpes (UGA).

[Zayas-Cabán et al., 2017]
Zayas-Cabán G, Jasin S, Wang G (2017) An asymptotically optimal heuristic
for general non-stationary finite-horizon restless multi-armed multi-action
bandits. Ross: Technology & Operations (Topic) .

[Zhang  Frazier, 2021]
Zhang X, Frazier PI (2021) Restless bandits with many arms: Beating the central
limit theorem. arXiv preprint arXiv:2107.11911 .

[Zipkin, 2000]
Zipkin PH (2000) Foundations of inventory management.

Certainty Equivalence Control-Based Heuristics in Multi-Stage Convex Stochastic Optimization Problems (Online Appendix)

*Outline Appendices appendix:proof-well-definedEC.1 and appendix:proof-of-main-resultsEC.2 contain the technical proofs of our theoretical results. In Appendix appendix:policy-mappingEC.3, we provide a unifying view on the various regularity conditions seen as sufficient for specific upper bounds. Appendix appendix:policy-classes-extensionEC.4 extends our proposed heuristic policies to several broader policy classes, followed by a discussion on the issue of the multiplicative exponentially growing constant.

§ EXISTENCE OF OPTIMAL SOLUTION

In this appendix, we provide a proof for Proposition <ref>, thereby validating the well-posedness of the optimization problem outlined in (<ref>). A key aspect within the proof process involves establishing the continuity of a maximization function with respect to a variable that is also part of the constraint set. To ensure this continuity, we leverage a technical result from Proposition 4.4 in [Bonnans  Shapiro, 2013], which affirms this continuity provided certain constraint qualification is met on the constrained sets.

[Proposition 1]
Under Assumptions <ref>-<ref>, the optimization problem (<ref>) is well-defined and an optimal solution exists. Moreover, the mapping $\bx \mapsto V_{\mathrm{opt}} (\bx,T)$ is a continuous function of $\bx$.

Proof of Proposition <ref>
We use dynamic programming and proceed backward in time. Starting at the last time-step $T$, define for $(\bx,\bw)$ the optimal value function

\begin{equation*}
\hat{V}_{\mathrm{opt}} (\bx,1,\bw) \defeq \max_{\bu \in \calU_T (\bx,\bw)} R_T (\bx, \bw, \bu)
\end{equation*}

By Weierstrass extreme value theorem, $\hat{V}_{\mathrm{opt}} (\bx,1,\bw)$ is well-defined and attained, since it is the maximum value of a continuous function $ \bu \mapsto R_T(\bx, \bw, \bu)$ over a non-empty compact set $\calU_T (\bx,\bw)$, where the compactness follows from our model assumptions.

We next show that $\bx \mapsto \hat{V}_{\mathrm{opt}} (\bx,1,\bw)$ is continuous as a function of $\bx$, by using <cit.>. Note that this proposition is established under a more general setting of infinite dimensional Banach space, and for our finite dimensional case considered here, the only non-trivial condition that needs verification is Robinson's Constraint Qualification (Condition (iv) of this Proposition 4.4). In finite dimension, Robinson's CQ reduces to MFCQ (Mangasarian-Fromovitz Constraint Qualification, see also Appendix appendix:policy-mappingEC.3 for more details). By <cit.>, the MFCQ holds at any point $\bu \in \calU_T(\bx,\bw)$ if and only if the Slater CQ holds for $\calU_T(\bx,\bw)$. By the first item of Assumption <ref>, this is true for any feasible set $\calU_T(\bx,\bw)$ parameterized by $(\bx,\bw)$. We conclude that $\bx \mapsto \hat{V}_{\mathrm{opt}} (\bx,1,\bw)$ is indeed a continuous function of $\bx$. Moreover, since $\bx$ takes values in a compact set, we deduce that $\bx \mapsto \hat{V}_{\mathrm{opt}} (\bx,1,\bw)$ is uniformly continuous.

By definition, we have

\begin{equation*}
V_{\mathrm{opt}} (\bx,1) = \int \hat{V}_{\mathrm{opt}} (\bx,1,\bw) \cdot f(\bw) d \bw
\end{equation*}

Since $\bx \mapsto \hat{V}_{\mathrm{opt}} (\bx,1,\bw)$ is uniformly continuous for any realisation of $\bw$, we deduce that $\bx \mapsto V_{\mathrm{opt}} (\bx,1)$ is also a continuous function of $\bx$.

Next for time-step $T-1$, the dynamic programming equation writes

\begin{align*}
\hat{V}_{\mathrm{opt}} (\bx,2,\bw) & = \max_{\bu \in \calU_{T-1} (\bx,\bw)} R_{T-1} (\bx, \bw, \bu) + \expect{V_{\mathrm{opt}} (\bX,1) \ \Big | \ \bx, \bw, \bu} \\
& \ = \max_{\bu \in \calU_{T-1} (\bx,\bw)} R_{T-1} (\bx, \bw, \bu) + \int_{\bx' \in \R^{n_x}} V_{\mathrm{opt}} (\phi(\bx,\bw,\bu) + \bx' ,1) \ d \nu(\bx' \mid \bx, \bw, \bu)
\end{align*}

We argue that $\bu \mapsto R_{T-1} (\bx, \bw, \bu) + \int_{\bx' \in \R^{n_x}} V_{\mathrm{opt}} (\phi(\bx,\bw,\bu) + \bx' ,1) \ d \nu(\bx' \mid \bx, \bw, \bu)$ is a continuous function of $\bu$, by relying on the second item of Assumption <ref> for the continuity of the mapping $\bu \mapsto \nu(\cdot \mid \bx,\bw,\bu)$, as well as the fact that $\bu \mapsto R_{T-1} (\bx, \bw, \bu)$ and $\bu \mapsto V_{\mathrm{opt}} (\phi(\bx,\bw,\bu) + \bx' ,1)$ are all continuous, the later follows from our previous analysis for time-step $T$. Consequently, again by Weierstrass extreme value theorem, $\hat{V}_{\mathrm{opt}} (\bx,2,\bw)$ is well-defined and attained for every $(\bx,\bw)$. We show in a similar way as for time-step $T$ that, $\bx \mapsto R_{T-1} (\bx, \bw, \bu) + \int_{\bx' \in \R^{n_x}} V_{\mathrm{opt}} (\phi(\bx,\bw,\bu) + \bx' ,1) \ d \nu(\bx' \mid \bx, \bw, \bu)$ is a continuous function of $\bx$. Again by the first item of Assumption <ref> and <cit.>, we deduce that $\bx \mapsto \hat{V}_{\mathrm{opt}} (\bx,2,\bw)$ is a (uniformly) continuous function of $\bx$ for any realisation of $\bw$. Consequently $\bx \mapsto V_{\mathrm{opt}} (\bx,2) \defeq \int \hat{V}_{\mathrm{opt}} (\bx,2,\bw) \cdot f(\bw) d \bw$ is also continuous.

More generally, suppose that for time-step $t$ the optimization problem $V_{\mathrm{opt}} (\bx,t)$ of (<ref>) for time-span $[T-t+1,T]$ is well-defined and $\bx \mapsto V_{\mathrm{opt}} (\bx,t)$ is a continuous function of $\bx$. We remind the reader that since we are using backward induction, the "$t$" in $V_{\mathrm{opt}} (\bx,t)$ is indexing the time steps remaining until the end of horizon. We write the dynamic programming equation

\begin{equation*}
\hat{V}_{\mathrm{opt}} (\bx,t+1,\bw) = \max_{\bu \in \calU_{T-t} (\bx,\bw)} R_{T-t} (\bx, \bw, \bu) + \int_{\bx' \in \R^{n_x}} V_{\mathrm{opt}} (\phi(\bx,\bw,\bu) + \bx', t) \ d \nu(\bx' \mid \bx, \bw, \bu)
\end{equation*}

We use the induction hypothesis and item 2 of Assumption <ref> to show that $\bu \mapsto R_{T-t} (\bx, \bw, \bu) + \int_{\bx' \in \R^{n_x}} V_{\mathrm{opt}} (\phi(\bx,\bw,\bu) + \bx', t) \ d \nu(\bx' \mid \bx, \bw, \bu)$ is continuous, so by Weierstrass extreme value theorem $\hat{V}_{\mathrm{opt}} (\bx,t+1,\bw)$ is wel-defined and attained. We use again the induction hypothesis and item 2 of Assumption <ref> to show that $\bx \mapsto R_{T-t} (\bx, \bw, \bu) + \int_{\bx' \in \R^{n_x}} V_{\mathrm{opt}} (\phi(\bx,\bw,\bu) + \bx', t) \ d \nu(\bx' \mid \bx, \bw, \bu)$ is continuous, and combine with item 1 of Assumption <ref> to show that $\bx \mapsto \hat{V}_{\mathrm{opt}} (\bx,t+1,\bw)$ is (uniformly) continuous for any realisation of $\bw$. Hence $\bx \mapsto V_{\mathrm{opt}} (\bx,t+1) \defeq \int \hat{V}_{\mathrm{opt}} (\bx,t+1,\bw) \cdot f(\bw) d \bw$ is continuous as well. This completes the induction step and concludes the proof.

To establish the continuity result presented in Proposition <ref>, we relied on the MFCQ. As we will elaborate in Appendix appendix:policy-mappingEC.3, stronger constraint qualifications are required to ensure stricter regularities, such as Lipschitz-continuity or $\calC^2$-smoothness.

§ PROOF OF THE PERFORMANCE BOUNDS

In this appendix we prove the performance bounds in Section <ref>. The following concentration inequality in vector form will be used in the sequence to bound the probability of the stochastic trajectory leaving outside an $\varepsilon$-neighbourhood.

Let $\bY$ be a random vector such that $\expect{\bY} = \by$, $\norme{\bY} < C$ and $\var{\bY} < \infty$. Then for $0 < \varepsilon < \var{\bY}/C$, we have

\begin{equation*}
\proba{\norme{\bY - \by} \ge \varepsilon} \le \exp \left( - \frac{\varepsilon^2}{8 \var{\bY}} + \frac{1}{4} \right)
\end{equation*}

For a proof of Lemma <ref> we refer to <cit.>.

§.§ The Update Policy
For convenience, we repeat the assumptions and theorems below for ease of discussion.

[Assumption I]
For all $1 \le t \le T$, there exists $\varepsilon_t > 0$ such that for all $(\bx,\bw) \in \calB((\bx^*(t),\barw), \varepsilon_t)$, the set $S^*_t(\bx,\bw)$ is single-valued. Moreover, the (single-valued) function $S^*_t(\bx,\bw)$ defined in (<ref>) is locally Lipschitz-continuous in $\calB((\bx^*(t),\barw), \varepsilon_t)$.

Combining the smoothness (hence the Lipschitz-continuity) of the reward functions $R_t(\cdot)$, Assumption <ref> implies that

* The function $\hat{V}_{\mathrm{rel}+} (\cdot, T-t, \cdot): \R^{n_x} \times \R^{n_w} \rightarrow \R$ defined in (<ref>) is Lipschitz-continuous with Lipschitz constant $K_t$ in the domain $\calB((\bx^*(t),\barw), \varepsilon_t)$:

\begin{equation} \label{eq:lipschitz-approximation-1}
\abs{\hat{V}_{\mathrm{rel}+} (\bx, T-t, \bw) - \hat{V}_{\mathrm{rel}+} (\bx^*(t), T-t, \barw)} \le K_t \cdot \norme{(\bx, \bw) - (\bx^*(t), \barw)}, \ \ \ \forall (\bx,\bw) \in \calB((\bx^*(t),\barw), \varepsilon_t)
\end{equation}

Since $V_{\mathrm{rel}-} (\bx,T+1-t) = \hat{V}_{\mathrm{rel}+} (\bx,T+1-t,\barw)$, this also implies that

* The function $V_{\mathrm{rel}-} (\cdot, T-t): \R^{n_x} \rightarrow \R$, defined in (<ref>) is Lipschitz-continuous with Lipschitz constant $L_t$ in the domain $\calB(\bx^*(t), \varepsilon_t)$:

\begin{equation}  \label{eq:lipschitz-approximation-2}
\abs{ V_{\mathrm{rel}-} (\bx,T+1-t) - V_{\mathrm{rel}-} (\bx^*(t),T+1-t) } \le L_t \cdot \norme{\bx - \bx^*(t)}, \ \ \ \forall \bx \in \calB(\bx^*(t), \varepsilon_t)
\end{equation}

[Theorem 1 (restated with explicit constants)]
Let $V_{\mathrm{opt}} (\bx,T)$ be the value of the stochastic optimization problem (<ref>) that satisfies Assumptions <ref>-<ref>, and let $V_{\mathrm{update}} (\bx,T)$ be the value of the update policy defined in Algorithm <ref>. Under the additional Assumption <ref>, there exists constants $C_1, \mathfrak{P}, \bar{V} > 0$ such that

\begin{equation*}
V_{\mathrm{opt}} (\bx,T) - V_{\mathrm{update}} (\bx,T) \le \mathfrak{P} C_1 +  (1 - \mathfrak{P}) \bar{V}
\end{equation*}

The constant $\mathfrak{P}$ is given in (<ref>) and converges to $1$ exponentially fast as both $\var{\bW}$ and $\var{\calE}$ converge to $0$. The constant $\bar{V}$ is a finite upper bound of $V_{\mathrm{opt}} (\bx,T)$.
The constant $C_1$ is

\begin{equation}  \label{eq:choice-of-constant}
C_1 \defeq \sum_{t=1}^{T} K_{t} \sqrt{\var{\bW}} + L_{t} \sqrt{\var{\calE}}
\end{equation}

* $K_t$ is the Lipschitz constant of the function $\hat{V}_{\mathrm{rel}+} (\cdot, T-t, \cdot): \R^{n_x} \times \R^{n_w} \rightarrow \R$ defined in (<ref>) for $(\bx,\bw) \in \calB((\bx^*(t),\barw), \varepsilon_t)$

* $L_t$ is the Lipschitz constant of the function $V_{\mathrm{rel}-} (\cdot, T-t): \R^{n_x} \rightarrow \R$ defined in (<ref>) for $\bx \in \calB(\bx^*(t), \varepsilon_t)$

* The distribution function of the random vector $\bW$ are given in Assumption <ref>

* $\var{\calE} < \infty$ is defined in (<ref>) of Assumption <ref>

For each time-step $1 \le t \le T$, we introduce the following value function

\begin{equation} \label{eq:rel-bar-t}
V_{\mathrm{rel}+} (\bx(t),T+1-t) \defeq \int \hat{V}_{\mathrm{rel}+} (\bx(t), T+1-t, \bw) \cdot f(\bw) d \bw
\end{equation}

where recall that $f(\bw)$ is the distribution function of the random vector $\bW$. An interpretation of $V_{\mathrm{rel}+} (\bx(t),T+1-t)$ is that it is the expected value of the deterministic optimal control at time-step $t$, with the system being in $\bx(t)$, plus the additional hindsight of the exact value of $\bw = \bW(t)$. In contrast, $V_{\mathrm{rel}-} (\bx(t),T+1-t)$ from (<ref>) can be seen as the expected value of the deterministic optimal control at time-step $t$, with the system being in $\bx(t)$, and zero knowledge of the realization of $\bW(t)$. So we have the bound $V_{\mathrm{rel}-} (\bx(t),T+1-t) \le V_{\mathrm{rel}+} (\bx(t),T+1-t)$. Their difference can be interpreted as the value of perfect information of knowing $\bW$, see [Avriel  Williams, 1970, Huang et al., 1977]. We are now ready to prove the theorem.

Proof of Theorem <ref>
We divide the proof into several steps.

Step One: Bounding the Probability of Leaving the $\varepsilon_t$-Neighbourhood
For $1 \le t \le T$, denote by $J_t$ the Lipschitz constant of the policy mapping $S^*_t$ inside $\calB((\bx^*(t),\barw), \varepsilon_t)$ claimed in Assumption <ref>, and write $J \defeq \max_t J_t$. We have

\begin{equation*}
\bx^*(t+1) = \phi(\bx^*(t), \barw, S^*_t(\bx^*(t), \barw))
\end{equation*}

\begin{equation*}
\bX(t+1) = \phi(\bX(t), \bW(t), S^*_t(\bX(t), \bW(t))) + \calE(t)
\end{equation*}

where we have abbreviated $\calE(t) \defeq \calE(\bX(t), \bW(t), S^*_t(\bX(t), \bW(t)))$. Since $\phi(\cdot)$ is an affine function, we denote by $c_{\phi}$ its Lipschitz constant. Making the difference of the above two equations and re-arranging terms, we obtain

\begin{equation*}
\norme{\bX(t+1) - \bx^*(t+1)} \le c_{\phi} (J+1) \norme{\bX(t) - \bx^*(t)} + c_{\phi} (J+1) \norme{\bW(t) - \barw} + \norme{\calE(t)}
\end{equation*}

Denote by $a \defeq c_{\phi} (J+1)$, $\varepsilon \defeq \min_t \varepsilon_t/2$, and

\begin{equation*}
b \defeq \begin{cases}
\varepsilon/T, & \mbox{if } a = 1 \\
\varepsilon(a-1)/a^T, & \mbox{if } a > 1 \\
\varepsilon (1 - a), & \mbox{if } 0 < a < 1
\end{cases}
\end{equation*}

An elementary calculation shows that

\begin{equation*}
\mbox{$a \norme{\bW(t) - \barw} + \norme{\calE(t)} \le b$ holds for all $t$} \Rightarrow \norme{\bX(t) - \bx^*(t)} \le \varepsilon_t/2 \mbox{ holds for all } t
\end{equation*}

Hence we deduce that

\begin{align}
& \norme{\calE(t)} \le b/2 \mbox{ and } \norme{\bW(t) - \barw} \le \min \{ b/(2a), \varepsilon \} \mbox{ hold for all time-step $t$}  \nonumber \\
\Rightarrow & \qquad \mbox{the stochastic trajectory $\bX(t)$ together with all the realizations of $\bW(t)$ during the update} \nonumber \\
& \qquad \mbox{policy remains inside the $\varepsilon_t$-neighbourhood required by Assumption \ref{assumption:lipschitz-continuity} for all time-step $t$}  \label{eq:event-for-time-t}
\end{align}

By Lemma <ref>, the event described on the left hand side of (<ref>) occurs with probability at least

\begin{equation} \label{eq:event-proba-formula}
\mathfrak{P} \defeq \left[1 - \exp \left(- \frac{(b/2)^2}{8 \var{\calE}} + \frac{1}{4} \right) \right ]^T  \left[1 - \exp \left(- \frac{(\min \{ b/(2a), \varepsilon \})^2}{8 \var{\bW}} + \frac{1}{4} \right) \right ]^T
\end{equation}

Step Two: Computation Inside the $\varepsilon_t$-Neighbourhood
In this step we suppose that the event described in (<ref>) occurs. Fix $1 \le t \le T$ and suppose that arriving at time-step $t$, the system configuration is in $\bX(t) = \bx(t)$. We can bound the difference between $V_{\mathrm{rel}-} (\bx(t),T+1-t)$ and $V_{\mathrm{rel}+} (\bx(t),T+1-t)$ by

\begin{align} \label{eq:difference-of-prio-and-post}
& \abs{V_{\mathrm{rel}-} (\bx(t),T+1-t) - V_{\mathrm{rel}+} (\bx(t),T+1-t)} \\
\le &   \int \abs{ \hat{V}_{\mathrm{rel}+} (\bx(t), T+1-t, \expect{\bW}) - \hat{V}_{\mathrm{rel}+} (\bx(t), T+1-t, \bw) } \cdot f(\bw) d \bw \nonumber \\
\le & \ K_{t-1} \int \norme{\expect{\bW} - \bw} \cdot f(\bw) d \bw \le K_{t-1} \sqrt{\var{\bW}}  \nonumber
\end{align}

where in the last line we have applied (<ref>) from Assumption <ref> with the implication from event (<ref>). Next observe that according to the definition of the update policy, it is planned by using the full knowledge of $\bw = \bW(t)$, so we can write

\begin{equation} \label{eq:update-bar-t}
V_{\mathrm{update}} (\bx(t),T+1-t) = \int \hat{V}_{\mathrm{update}} (\bx(t), T+1-t, \bw) \cdot f(\bw) d \bw
\end{equation}

\begin{align}
& \hat{V}_{\mathrm{update}} (\bx(t), T+1-t, \bw) \defeq R_{t} \left( \bx(t), \bw, \bu^*_{\bx,T+1-t,\bw}(t) \right)  \nonumber \\
& + \expect{V_{\mathrm{update}} (\bX(t+1),T-t) \ \Big| \ \bX(t) = \bx(t), \bW(t) = \bw, \bU(t) = \bu^*_{\bx(t),T+1-t,\bw}(t)} \label{eq:principle-of-optimality-1}
\end{align}

In the above formula, $R_{t} \left( \bx(t), \bw, \bu^*_{\bx,T+1-t,\bw}(t) \right)$ is the instantaneous reward gained at time-step $t$, with the system being in configuration $\bx(t)$ and the control $\bu^*_{\bx,T+1-t,\bw}(t)$ obtained by solving (<ref>) is applied. The system then evolves to the configuration $\bX(t+1)$ at time-step $t+1$, and the second term accounts for the expected reward onwards.

On the other hand, by Bellman's principle of optimality, we claim that

\begin{equation}  \label{eq:principle-of-optimality}
\hat{V}_{\mathrm{rel}+} (\bx(t),T+1-t,\bw) = R_{t} \big( \bx(t), \bw, \bu^*_{\bx(t),T+1-t,\bw}(t) \big) +  V_{\mathrm{rel}-} (\phi \big(\bx(t), \bw, \bu^*_{\bx(t),T+1-t,\bw}(t) \big),T-t)
\end{equation}

Indeed, $\bu^*_{\bx(t),T+1-t,\bw}(t)$ is an optimal control in the first-step from solving (<ref>) for $\hat{V}_{\mathrm{rel}+} (\bx(t),T+1-t,\bw)$, and $\hat{\bx}^*(t+1) \defeq \phi \big(\bx(t), \bw, \bu^*_{\bx(t),T+1-t,\bw}(t) \big)$ is an optimal trajectory in the second-step. Since from the second-step onward, the optimal control and the optimal trajectory for $\hat{V}_{\mathrm{rel}+} (\bx(t),T+1-t,\bw)$ coincide with the ones for $V_{\mathrm{rel}-} (\hat{\bx}^*(t+1),T-t)$ in (<ref>), we deduce (<ref>).

We remark that the reward at time-step $t$ in $\hat{V}_{\mathrm{update}} (\bx(t), T+1-t, \bw)$ and $\hat{V}_{\mathrm{rel}+} (\bx(t),T+1-t,\bw)$ are the same. Hence by taking the difference of (<ref>) and (<ref>), and using (<ref>) and (<ref>), we obtain

\begin{align}  \label{eq:difference}
& V_{\mathrm{update}} (\bx(t),T+1-t) - V_{\mathrm{rel}+} (\bx(t),T+1-t) \nonumber  \\
=  & \int \expect{V_{\mathrm{update}} (\bX(t+1),T-t) \ \Big| \ \bX(t) = \bx(t), \bW(t) = \bw, \bU(t) = \bu^*_{\bx(t),T+1-t,\bw}(t)} \cdot f(\bw) d \bw \nonumber \\
& \qquad - \int V_{\mathrm{rel}-} (\phi \big(\bx(t), \bw, \bu^*_{\bx(t),T+1-t,\bw}(t) \big),T-t) \cdot f(\bw) d \bw
\end{align}

We write out the following decomposition, using the abbreviation $\bx^*(t+1) \defeq \phi \big(\bx(t), \bw, \bu^*_{\bx(t),T+1-t,\bw}(t) \big)$ and $\bu^*(t) \defeq \bu^*_{\bx(t),T+1-t,\bw}(t)$:

\begin{align} \label{eq:decomposition-of-difference}
& \expect{V_{\mathrm{update}} (\bX(t+1),T-t) \ \Big| \ \bX(t) = \bx(t), \bW(t) = \bw, \bU(t) = \bu^*(t)} - V_{\mathrm{rel}-} \big(\bx^*(t+1),T-t \big) \nonumber   \\
= &  \ \underbrace{\expect{ V_{\mathrm{update}} (\bX(t+1),T-t) - V_{\mathrm{rel}+} (\bX(t+1),T-t) \ \Big| \ \bX(t) = \bx(t), \bW(t) = \bw, \bU(t) = \bu^*(t)}}_{\mytag{Term A}{term:A}} \nonumber \\
& + \underbrace{\expect{ V_{\mathrm{rel}+} (\bX(t+1),T-t) - V_{\mathrm{rel}-} (\bX(t+1),T-t) \ \Big| \ \bX(t) = \bx(t), \bW(t) = \bw, \bU(t) = \bu^*(t)}}_{\mytag{Term B}{term:B}}  \\   \nonumber
& + \underbrace{\expect{ V_{\mathrm{rel}-} (\bX(t+1),T-t) - V_{\mathrm{rel}-} \big(\bx^*(t+1),T-t \big) \ \Big| \ \bX(t) = \bx(t), \bW(t) = \bw, \bU(t) = \bu^*(t) }}_{\mytag{Term C}{term:C}}   \nonumber
\end{align}

We now analyse the three terms in (<ref>). Denote by

\begin{equation*}
Z(t,\bX(t)) \defeq \expect{V_{\mathrm{update}} (\bX(t),T+1-t) - V_{\mathrm{rel}+} (\bX(t),T+1-t) \ \Big| \ \bX(t)}
\end{equation*}

which is interpreted as follows: conditional on the value of $\bX(t)$, the difference between the expected value of the update policy performance on system state $\bX(t)$ for the time-span $[t , T]$ with the upper bound $V_{\mathrm{rel}+} (\bX(t),T+1-t)$. Applying $\expect{\ \cdot \ \big| \ \bX(t)}$ in (<ref>) and use (<ref>), we obtain

\begin{equation} \label{eq:decomposition-of-terms}
Z(t,\bX(t)) = \expect{\int (\mbox{\ref{term:A}} + \mbox{\ref{term:B}} + \mbox{\ref{term:C}} ) \cdot f(\bw) d \bw \ \Big| \ \bX(t)}
\end{equation}

For $\mbox{\ref{term:A}}$, we have

\begin{align}  \label{eq:analysis-of-term-A}
& \expect{\int (\mbox{\ref{term:A}} ) \cdot f(\bw) d \bw \ \Big| \ \bX(t)}  \nonumber \\
= & \ \expect{\int \expect{ V_{\mathrm{update}} (\bX(t+1),T-t) - V_{\mathrm{rel}+} (\bX(t+1),T-t) \ \Big| \ \bW(t) = \bw, \bU(t) = \bu^*(t)}  \cdot f(\bw) d \bw \ \Big | \ \bX(t) }  \nonumber \\
= & \ \expectup{ Z(t+1,\bX(t+1)) \ \Big| \ \bX(t)}
\end{align}

where $\expectup{ Z(t+1,\bX(t+1)) \ \big| \ \bX(t)} $ means conditional on $\bX(t)$, the expected value of the quantity $Z(t+1,\bX(t+1))$ under the update policy.

For $\mbox{\ref{term:B}}$, by (<ref>) we have

\begin{equation*}
\norme{ \expect{\int (\mbox{\ref{term:B}} ) \cdot f(\bw) d \bw} \ \Big| \ \bX(t) } \le K_{t} \sqrt{\var{\bW}}
\end{equation*}

For $\mbox{\ref{term:C}}$, we have

\begin{align*}
& \norme{ \expect{\int (\mbox{\ref{term:C}} ) \cdot f(\bw) d \bw} \ \Big| \ \bX(t) } \\
\le & \ L_{t} \int \expect{\norme{\bX(t+1) - \bx^*(t+1)} \ \big| \ \bX(t) = \bx(t), \bW(t) = \bw, \bU(t) = \bu^*(t) } \cdot f(\bw) d \bw  \\
\le & \ L_{t} \sqrt{\var{\calE}}
\end{align*}

where in the first inequality we have applied (<ref>) from Assumption <ref> with the implication from event (<ref>), and the second inequality follows from Assumption <ref>.

Now if we write

\begin{equation} \label{eq:definition-of-small-z}
z(t) \defeq \expectup{Z(t,\bX(t)) \mid \bX(1) = \bx}
\end{equation}

which is the expected value under the update policy for the quantity $Z(t,\bX(t))$ conditional merely on the initial system state $\bX(1)$, we have

\begin{align*}
z(t) - z(t+1) &  = \expectup{Z(t,\bX(t)) - Z(t+1, \bX(t+1)) \ \big| \ \bX(1) = \bx } \\
& = \expectup{Z(t,\bX(t)) - \expectup{Z(t+1, \bX(t+1)) \mid  \bX(t) }  \ \big| \ \bX(1) = \bx }
\end{align*}

So from the above calculations we obtain

\begin{equation*}
\abs{z(t) - z(t+1)}  \le K_{t} \sqrt{\var{\bW}} + L_{t} \sqrt{\var{\calE}}
\end{equation*}

Step Three: Conclusion of the Proof

\begin{equation*}
z(1) = V_{\mathrm{update}} (\bx,T) - V_{\mathrm{rel}+} (\bx,T)
\end{equation*}

and $z(T+1) = 0$, we conclude that

\begin{align*}
& \abs{V_{\mathrm{update}} (\bx,T) - V_{\mathrm{opt}} (\bx,T)} \le \abs{V_{\mathrm{update}} (\bx,T) - V_{\mathrm{rel}+} (\bx,T)}  \\
= & \  \abs{z(1)} \le \sum_{t=1}^{T} \abs{z(t) - z(t+1)}  \\
\le & \ \mathfrak{P} \left( \sum_{t=1}^{T} K_{t} \sqrt{\var{\bW}} + L_{t} \sqrt{\var{\calE}} \right) + (1 - \mathfrak{P}) \bar{V}
\end{align*}

We next recall Assumption <ref>:

[Assumption II]
For all $1 \le t \le T$, there exists $\varepsilon_t > 0$ such that for all $(\bx,\bw) \in \calB((\bx^*(t),\barw), \varepsilon_t)$, the set $S^*_t(\bx,\bw)$ is single-valued. Moreover, the (single-valued) function $S^*_t(\bx,\bw)$ defined in (<ref>) is locally $\calC^2$-smooth in $\calB((\bx^*(t),\barw), \varepsilon_t)$.

Combining the smoothness of the reward functions $R_t(\cdot)$, Assumption <ref> implies that

* The function $\hat{V}_{\mathrm{rel}+} (\cdot, T-t, \cdot): \R^{n_x} \times \R^{n_w} \rightarrow \R$ defined in (<ref>) is $\calC^2$-smooth in the domain $\calB((\bx^*(t),\barw), \varepsilon_t)$, so that

\begin{align} \label{eq:second-order-approximation-1}
\hat{V}_{\mathrm{rel}+} (\bx, T-t, \bw) & = \hat{V}_{\mathrm{rel}+} (\bx^*(t), T-t, \barw) + \left[(\bx,\bw) - (\bx^*(t),\barw)\right] \cdot \nabla_{(\bx,\bw)}^\top \hat{V}_{\mathrm{rel}+} (\bx^*(t), T-t, \barw)  \nonumber \\
& \qquad + \frac{1}{2} \left[(\bx,\bw) - (\bx^*(t),\barw)\right] \cdot \nabla^2_{(\bx,\bw)} \hat{V}_{\mathrm{rel}+} (\bx', T-t, \bw') \cdot \left[(\bx,\bw) - (\bx^*(t),\barw)\right]^\top
\end{align}

for a certain $(\bx',\bw')$ along the line segment from $(\bx, \bw)$ to $(\bx^*(t),\barw)$ in $\R^{n_x} \times \R^{n_w}$. By our choice of norms, we have

\begin{align*}
& \abs{\left[(\bx,\bw) - (\bx^*(t),\barw)\right] \cdot \nabla^2_{(\bx,\bw)} \hat{V}_{\mathrm{rel}+} (\bx', T-t, \bw') \cdot \left[(\bx,\bw) - (\bx^*(t),\barw)\right]^\top} \\
& \qquad \le \norme{(\bx,\bw) - (\bx^*(t),\barw)}^2 \norme{\nabla^2_{(\bx,\bw)} \hat{V}_{\mathrm{rel}+} (\bx', T-t, \bw')}
\end{align*}

Denote by $\mathfrak{K}_t (\bx',\bw') \defeq \norme{\nabla^2_{(\bx,\bw)} \hat{V}_{\mathrm{rel}+} (\bx', T-t, \bw')}$. We have

\begin{equation*}
\sup_{(\bx', \bw') \in \calB((\bx^*(t),\barw), \varepsilon_t)} \mathfrak{K}_t (\bx',\bw') \defeq \norminf{\mathfrak{K}_t} < \infty
\end{equation*}

* The function $V_{\mathrm{rel}-} (\cdot, T-t): \R^{n_x} \rightarrow \R$, defined in (<ref>) is $\calC^2$-smooth in the domain $\calB(\bx^*(t), \varepsilon_t)$, so that

\begin{align} \label{eq:second-order-approximation-2}
V_{\mathrm{rel}-} (\bx, T-t) & = V_{\mathrm{rel}-} (\bx^*(t), T-t) + (\bx - \bx^*(t)) \cdot \nabla_{\bx}^\top V_{\mathrm{rel}-} (\bx^*(t), T-t)  \nonumber \\
& \qquad + \frac{1}{2} (\bx - \bx^*(t)) \cdot \nabla^2_{\bx} V_{\mathrm{rel}-} (\bx', T-t) \cdot (\bx - \bx^*(t))^\top
\end{align}

for a certain $\bx'$ along the line segment from $\bx$ to $\bx^*(t)$ in $\R^{n_x}$. We have

\begin{equation*}
\sup_{\bx' \in \calB(\bx^*(t), \varepsilon_t)} \norme{\nabla^2_{\bx} V_{\mathrm{rel}-} (\bx', T-t)} \defeq \norminf{\mathfrak{L}_t} < \infty
\end{equation*}

[Theorem 2 (restated with explicit constants)]
Let $V_{\mathrm{opt}} (\bx,T)$ be the value of the stochastic optimization problem (<ref>) that satisfies Assumptions <ref>-<ref>, and let $V_{\mathrm{update}} (\bx,T)$ be the value of the update policy defined in Algorithm <ref>. Under the additional Assumption <ref>, there exists constants $C_2, \mathfrak{P}, \bar{V} > 0$ such that

\begin{equation*}
V_{\mathrm{opt}} (\bx,T) - V_{\mathrm{update}} (\bx,T) \le \mathfrak{P} C_2 +  (1 - \mathfrak{P}) \bar{V}
\end{equation*}

The constant $\mathfrak{P}$ is given in (<ref>) and converges to $1$ exponentially fast as both $\var{\bW}$ and $\var{\calE}$ converge to $0$. The constant $\bar{V}$ is a finite upper bound of $V_{\mathrm{opt}} (\bx,T)$.
The constant $C_2$ is

\begin{equation}  \label{eq:choice-of-constant-refined}
C_2 \defeq \frac{1}{2} \sum_{t=1}^{T} \norminf{\mathfrak{K}_t} \var{\bW} + \norminf{\mathfrak{L}_t} \var{\calE}
\end{equation}

* $\norminf{\mathfrak{K}_t} = \sup_{(\bx', \bw') \in \calB((\bx^*(t),\barw), \varepsilon_t)} \norme{\nabla^2_{(\bx,\bw)} \hat{V}_{\mathrm{rel}+} (\bx', T-t, \bw')}$

* $\norminf{\mathfrak{L}_t} = \sup_{\bx' \in \calB(\bx^*(t), \varepsilon_t)} \norme{\nabla^2_{\bx} V_{\mathrm{rel}-} (\bx', T-t)}$

* The distribution function of the random vector $\bW$ are given in Assumption <ref>

* $\var{\calE} < \infty$ is defined in (<ref>) of Assumption <ref>

Proof of Theorem <ref>
The proof uses exactly the same lines of analysis as in Theorem <ref>, up to arriving at Equation (<ref>). From there, using (<ref>), instead of (<ref>), we have

\begin{align} \label{eq:difference-of-prio-and-post-refine}
& V_{\mathrm{rel}-} (\bx(t),T+1-t) - V_{\mathrm{rel}+} (\bx(t),T+1-t)  \nonumber \\
= &   \int  \left( \hat{V}_{\mathrm{rel}+} (\bx(t), T+1-t, \expect{\bW}) - \hat{V}_{\mathrm{rel}+} (\bx(t), T+1-t, \bw) \right) \cdot f(\bw) d \bw \nonumber \\
= & \ \int (\expect{\bW}) - \bw) \cdot \nabla_{\bw}^\top \hat{V}_{\mathrm{rel}+} (\bx, T-t, \expect{\bW})) \cdot f(\bw) d \bw \nonumber \\
& \qquad + \frac{1}{2} \int \left( (\expect{\bW}) - \bw) \cdot \nabla^2_{\bw} \hat{V}_{\mathrm{rel}+} (\bx, T-t, \bw') \cdot (\expect{\bW}) - \bw)^\top  \right) \cdot f(\bw) d \bw \nonumber \\
= & \frac{1}{2} \int \left( (\expect{\bW}) - \bw) \cdot \nabla^2_{\bw} \hat{V}_{\mathrm{rel}+} (\bx, T-t, \bw') \cdot (\expect{\bW}) - \bw)^\top \right) \cdot f(\bw) d \bw
\end{align}

where in the above $\bw' = \bw'(\bw,\bw^*)$ depends on $\bw$ and $\bw^*$. The crucial part in (<ref>) is that the first-order term cancels out, and only the second-order term is left.

Similarly, using (<ref>) and the abbreviation $\bx^*(t+1) \defeq \phi \big(\bx(t), \bw, \bu^*_{\bx(t),T+1-t,\bw}(t) \big)$, $\calE(t) \defeq \bX(t+1) - \bx^*(t+1)$, our calculation concerning <ref> becomes

\begin{align*}
& V_{\mathrm{rel}-} (\bX(t+1),T-t) - V_{\mathrm{rel}-} \big(\bx^*(t+1),T-t \big)   \\
= & \   \calE(t) \cdot \nabla_{\bx}^\top V_{\mathrm{rel}-} (\bx^*(t+1), T-t) + \frac{1}{2} \calE(t) \cdot \nabla^2_{\bx} V_{\mathrm{rel}-} (\bx', T-t) \cdot \calE(t)^\top
\end{align*}

where in the above $\bx' = \bx'(\bX(t+1), \bx^*(t+1))$ depends on $\bX(t+1)$ and $\bx^*(t+1)$. Consequently

\begin{equation}  \label{eq:analysis-of-term-C}
\mbox{\ref{term:C}} = \frac{1}{2} \expectt{\calE(t) \cdot \nabla^2_{\bx} V_{\mathrm{rel}-} (\bx', T-t) \cdot \calE(t)^\top \ \big| \ \bX(t) = \bx(t), \bW(t) = \bw, \bU(t) = \bu^*_{\bx(t),T+1-t,\bw}(t)}
\end{equation}

Again the crucial part in (<ref>) is that the first-order term cancels out after taking expectation, and we are left with only the second-order term. Combining (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and using Assumption <ref>, we finally obtain

\begin{equation*}
\abs{z(t) - z(t+1)} \le \frac{1}{2} \norminf{\mathfrak{K}_t} \var{\bW} + \frac{1}{2} \norminf{\mathfrak{L}_t} \var{\calE}
\end{equation*}

The rest of the analysis are similar to Theorem <ref>.

§.§ The Projection Policy
For convenience, we repeat the assumptions and theorems below for ease of discussion.

[Assumption III]
For all $1 \le t \le T$, there exists $\varepsilon_t > 0$ such that $\proj_t(\bx,\bw)$ as defined in (<ref>) is locally Lipschitz-continuous in $\calB((\bx^*(t),\barw), \varepsilon_t)$.

[Theorem 3 (restated with explicit constants)]
Let $V_{\mathrm{opt}} (\bx,T)$ be the value of the stochastic optimization problem (<ref>) that satisfies Assumptions <ref>-<ref>, and let $V_{\mathrm{proj}} (\bx,T)$ be the value of the projection policy defined in Algorithm <ref>. Under the additional Assumption <ref>, there exists constants $C_3, \mathfrak{P}', \bar{V} > 0$ such that

\begin{equation*}
V_{\mathrm{opt}} (\bx,T) - V_{\mathrm{proj}} (\bx,T) \le \mathfrak{P}' C_3 +  (1 - \mathfrak{P}') \bar{V}
\end{equation*}

The constant $\mathfrak{P}'$ is given in (<ref>) and converges to $1$ exponentially fast as both $\var{\bW}$ and $\var{\calE}$ converge to $0$. The constant $\bar{V}$ is a finite upper bound of $V_{\mathrm{opt}} (\bx,T)$. The constant $C_3$ is

\begin{equation}  \label{eq:choice-of-constant-projection-Lipschitz}
C_3 \defeq \sum_{t=1}^{T} \left( c_{\bW}(t) \cdot \sqrt{\var{\bW}} + c_{\calE}(t) \cdot \sqrt{\var{\calE}} \right)
\end{equation}

\begin{equation*}
c_{\bW}(t) \defeq c_R (\cproj + 1) + \frac{c_R c_{\phi} (\cproj + 1)^2 ( c_{\phi}^{t-1} (\cproj + 1)^{t-1} - 1)}{ c_{\phi} (\cproj + 1) - 1}
\end{equation*}

\begin{equation*}
c_{\calE}(t) \defeq \frac{c_R (\cproj + 1) ( c_{\phi}^{t-1} (\cproj + 1)^{t-1} - 1)}{ c_{\phi} (\cproj + 1) - 1}
\end{equation*}

* $c_{\phi}$ being the Lipschitz constant of $\phi(\cdot)$

* $c_R \defeq \max_t c_{R_t}$ and $c_{R_t}$ being the Lipschitz constant of the reward function $R_t(\cdot)$

* $c_{\mathrm{proj}} \defeq \max_t c_{\mathrm{proj},t}$ with $c_{\mathrm{proj},t}$ being the Lipschitz constant of the projection mapping $\proj(\bx,\bw,t)$ inside $\calB((\bx^*(t),\barw), \varepsilon_t)$ claimed in Assumption <ref>

Proof of Theorem <ref>
The first step of the proof is very similar to the one of Theorem <ref>. Namely, write

\begin{equation} \label{eq:one-step-deterministic}
\bx^*(t+1) = \phi(\bx^*(t), \barw, \proj_t(\bx^*(t), \barw))
\end{equation}

\begin{equation}  \label{eq:one-step-stochastic}
\bX(t+1) = \phi(\bX(t), \bW(t), \proj_t(\bX(t), \bW(t)) ) + \calE(t)
\end{equation}

where we have abbreviated $\calE(t) \defeq \calE(\bX(t), \bW(t), \proj_t(\bX(t), \bW(t)) )$. Taking the difference we obtain

\begin{equation} \label{eq:x-relation}
\norme{\bX(t+1) - \bx^*(t+1)} \le c_{\phi} (c_{\mathrm{proj}}+1) \norme{\bX(t) - \bx^*(t)} + c_{\phi} (c_{\mathrm{proj}}+1) \norme{\bW(t) - \barw} + \norme{\calE(t)}
\end{equation}

where $c_{\phi}$ is the Lipschitz constant of $\phi(\cdot)$, and $c_{\mathrm{proj}} \defeq \max_t c_{\mathrm{proj},t}$, where for $1 \le t \le T$, $c_{\mathrm{proj},t}$ is the Lipschitz constant of the projection mapping $\proj(\bx,\bw,t)$ inside $\calB((\bx^*(t),\barw), \varepsilon_t)$ claimed in Assumption <ref>.

Denote by $\alpha \defeq c_{\phi} (c_{\mathrm{proj}} + 1)$, $\varepsilon \defeq \min_t \varepsilon_t/2$, and

\begin{equation*}
\beta \defeq \begin{cases}
\varepsilon/T, & \mbox{if } \alpha = 1 \\
\varepsilon(\alpha-1)/\alpha^T, & \mbox{if } \alpha > 1 \\
\varepsilon (1 - \alpha), & \mbox{if } 0 < \alpha < 1
\end{cases}
\end{equation*}

An elementary calculation shows that

\begin{equation*}
\mbox{$\alpha \norme{\bW(t) - \barw} + \norme{\calE(t)} \le \beta$ holds for all $t$} \Rightarrow \norme{\bX(t) - \bx^*(t)} \le \varepsilon_t/2 \mbox{ holds for all } t
\end{equation*}

Hence we deduce that

\begin{align}
& \norme{\calE(t)} \le \beta/2 \mbox{ and } \norme{\bW(t) - \barw} \le \min \{ \beta/(2 \alpha), \varepsilon \} \mbox{ hold for all time-step $t$}  \nonumber \\
\Rightarrow & \qquad \mbox{the stochastic trajectory $\bX(t)$ together with all the realizations of $\bW(t)$ during the projection} \nonumber \\
& \qquad \mbox{policy remains inside the $\varepsilon_t$-neighbourhood required by Assumption \ref{assumption:lipschitz-continuity-projection} for all time-step $t$}  \label{eq:event-for-time-t-1}
\end{align}

By Lemma <ref>, the event described on the left hand side of (<ref>) occurs with probability at least

\begin{equation} \label{eq:event-proba-formula-1}
\mathfrak{P}' \defeq \left[1 - \exp \left(- \frac{(\beta/2)^2}{8 \var{\calE}} + \frac{1}{4} \right) \right ]^T  \left[1 - \exp \left(- \frac{(\min \{ \beta/(2\alpha), \varepsilon \})^2}{8 \var{\bW}} + \frac{1}{4} \right) \right ]^T
\end{equation}

We subsequently assume that the right hand side of (<ref>) holds. Taking expectation in (<ref>) we obtain

\begin{equation}  \label{eq:illustration-3}
\expectproj{\norme{\bX(t+1) - \bx^*(t+1)}} \le \sqrt{\var{\calE}} + c_{\phi} (\cproj + 1) \sqrt{\var{\bW}} + c_{\phi} (\cproj + 1) \expectproj{\norme{\bX(t) - \bx^*(t)}}
\end{equation}

where $\expectproj{\cdot} = \expectproj{ \ \cdot \mid \bX(1) = \bx}$ is the expectation taken under the projection policy. Using the abbreviation $\mathfrak{h} \defeq \sqrt{\var{\calE}} + c_{\phi} (\cproj + 1) \sqrt{\var{\bW}}$, an elementary calculation implies that for all $t$:

\begin{equation*}
\expectproj{\norme{\bX(t) - \bx^*(t)}} \le \mathfrak{h} \cdot \frac{( c_{\phi}^{t-1} (\cproj + 1)^{t-1} - 1)}{ c_{\phi} (\cproj + 1) - 1}
\end{equation*}

where the right hand side is interpreted as $\mathfrak{h}(t-1)$ if $ c_{\phi} (\cproj + 1) = 1$.

By definition, we have

\begin{equation*}
V_{\mathrm{rel}-} (\bx,T) = \sum_{t=1}^{T} R_t (\bx^*(t),\barw,\bu^*(t))
\end{equation*}

\begin{equation*}
V_{\mathrm{proj}} (\bx,T) =  \expectproj{\sum_{t=1}^{T} R_t \left( \bX(t), \bW(t), \proj(\bX(t),\bW(t),t) \right)}
\end{equation*}

Let $c_{R_t}$ be the Lipschitz constant of the reward function $R_t(\cdot)$ and write $c_R \defeq \max_t c_{R_t}$. We deduce that

\begin{align}
& \abs{V_{\mathrm{rel}-} (\bx,T) - V_{\mathrm{proj}} (\bx,T)}  \nonumber \\
\le & \ \sum_{t=1}^{T} c_R (\cproj + 1) \expectproj{\norme{\bX(t) - \bx^*(t)}} + c_R (\cproj + 1) \sqrt{\var{\bW}} \label{eq:illustration-4} \\
\le & \ \sum_{t=1}^{T} c_R (\cproj + 1) \left[ \mathfrak{h} \cdot \frac{( c_{\phi}^{t-1} (\cproj + 1)^{t-1} - 1)}{ c_{\phi} (\cproj + 1) - 1}  +  \sqrt{\var{\bW}} \right] \nonumber
\end{align}

By rearranging terms we achieve the optimality gap bound claimed in the theorem.