Approximate Dynamic Programming: An Efficient Machine Learning Algorithm

Abstract

We propose an efficient machine learning algorithm for two-stage stochastic programs. This machine learning algorithm is termed as projected stochastic hybrid learning algorithm, and consists of stochastic sub-gradient and piecewise linear approximation methods. We use the stochastic sub-gradient and sample information to update the piecewise linear approximation on the objective function. Then we introduce a projection step, which implemented the sub-gradient methods, to jump out from a local optimum, so that we can achieve a global optimum. By the innovative projection step, we show the convergent property of the algorithm for general two-stage stochastic programs. Furthermore, for the network recourse problem, our algorithm can drop the projection steps, but still maintains the convergence property. Thus, if we properly construct the initial piecewise linear functions, the pure piecewise linear approximation method is convergent for general two-stage stochastic programs. The proposed approximate dynamic programming algorithm overcomes the high dimensional state variables using methods from machine learning, and its logic capture the critical ability of the network structure to anticipate the impact of decisions now on the future. The optimization framework, which is carefully calibrated against historical performance, make it possible to introduce changes in the decisions and capture the collective intelligence of the experienced decisions. Computational results indicate that the algorithm exhibits rapid convergence.

1. Introduction

Optimal learning addresses the challenge of how to collect information, as efficiently as possible, to make a decision in the present such that it minimizes the expectation of costs in the future with uncertainty. Collecting information is usually time consuming and expensive. For example, several large shippers, such as Amazon, Walmart, and IKEA, need to decide the quantity of products to ship from plants to warehouses to satisfy the retailers’ demand. The retailer usually makes their decisions before knowing the real demand. Then, after they know the retail demand, they optimize the shipping plans between retailers and warehouses. The aforementioned problems generally can be treated as the two-stage stochastic programs. The decisions that the retailer made now (Stage 1) will determines the state while solving the problem in future (Stage 2). Therefore, an optimal decision can be made now if we can compute the expected cost function (the recourse function) of Stage 2. In this chapter, we propose an efficient machine learning algorithm that can collect information very efficiently based on the knowledge gradient and solve the problem optimally. The main gap between MAT and proposed algorithm is that our proposed can collect the information based on knowledge gradient and overcomes the “curse of dimensionality”. Besides, it can transfer the problem into a polynomial solvable problem and has been proven convergent theoretically.

1.1 Motivation

Optimal learning is a rich filed that includes contributions from different communities. At the moment, this chapter focus on optimal learning in two-stage stochastic program, which is a practically important problem. Problems of this type arise in several areas in dynamic programming, in which the decision maker need to make temporal and spatial decisions before realizing events that will influence the decisions. For example, in empty container repositioning problems [1], shipping companies need to reposition empty containers before realizing the demand. In locomotive planning problems [2], railroads have to decide the schedule of trains in which locomotives are assigned before disruptions occur across the railway network. For relief distribution problems [3], humanitarian decision makers need to distribute emergency aid to disaster locations when the emergency aid materials are very scarce amidst great uncertainties. For job scheduling problems [4], the managers need to decide initial staffing levels and their working plans before the demand are realized. Most of the aforementioned applications are fully sequential problems, and they can be modeled as two-stage stochastic programming problems. Hence, the research of two-stage stochastic optimization in this chapter is very important. However, the main obstacle in most practical problem is that the expected cost function in Stage 2 is quite complex due to uncertainty. In this chapter, we propose a hybrid learning algorithm called projected stochastic hybrid learning algorithm (ProSHLA) to approximate the expected recourse function for two-stage stochastic programs. In order to demonstrate the efficiency of the algorithm, we also theoretically prove the convergence of the proposed algorithm mathematically.

In essence, ProSHLA is a hybrid of stochastic sub-gradient and piecewise linear approximation methods. The core of ProSHLA consists of a series of learning steps those provide information for updating the recourse function through a sequence of piecewise linear separable approximations, and a series of the projection steps those can guarantee convergence by implementing the stochastic sub-gradient method. The mathematical analysis and the computational results all demonstrates that when the initial piecewise linear approximation function is properly constructed for two-stage stochastic programs with network recourse, the learning algorithm can drop the projection steps without sacrificing convergence. Moreover, without the projection step, the learning algorithm only consists of a series of learning steps through a sequence of piecewise linear separable approximations, and can solve the practical complex problems very efficiently. Our innovative finding can help the practitioner and the scholar to understand the open problem that has puzzled them for decades: why does the piecewise linear approximation method can be efficient and convergent for stochastic programs with network recourse in practice. In this chapter, we provide the first theoretical support by our analytical results for the use of the piecewise linear approximation method in solving practical problems.

1.2 Literature review

In this chapter, we consider the two-stage stochastic programming problem as follows:

$$\mathit{min}{c}_0^{T}x+E_{\omega }\left[Q\left(x,\omega \right)\right]$$E1

s.t.,

$$\mathit{Ax}=b,$$

$$x\ge 0,$$

where $X\subset {\mathfrak{R}}^n$ denotes a convex compact set and the recourse function $Q\left(x,\omega \right)$ denotes the optimal value of the second stage problem:

$$Q\left(x,\omega \right)=min{c}_1^{T}y\left(\omega \right)$$E2

s.t.,

$$W\left(\omega \right)y=h\left(\omega \right)-T\left(\omega \right)x,$$

$$y\left(\omega \right)\ge 0\left(\omega \right).$$

In the above model, variables x and y denote the decision variables of stage 1 and 2 problems, respectively. $A$, $W\left(\omega \right)$ are constraint matrices, and parameters $c_0$ and $c_1$ denote the first and second stage vectors of cost coefficients, respectively.

Stochastic programming models and solution methods has been examined by many researchers. Comprehensive reviews and discussions were performed by Wallace and Ziemba [5]. The expected recourse function is extremely complex to evaluate except for a few special cases. There are various approximation methods those can be categorized into four groups. Let $\widehat{Q}\left(x\right)$ denote the approximate function. The first group includes scenario methods which use the sample average of $Q\left(x,{\omega }_i\right)$ for several samples, ${\omega }_1$,${\omega }_2$,…${\omega }_N$, to approximate the expected recourse function [6]. The approximation function is usually successively updated by the following function:

$$\widehat{Q}\left(x\right)=\frac{{\sum }_{i=1}^{N}Q\left(x,{\omega }_i\right)}{N}$$

Generally, the scenario method is very efficient, but it cannot guarantee to obtain the convergent solution.

The second group consists stochastic gradient techniques [7, 8], which updates solutions by using stochastic sub-gradients as directions. Usually, the approximate function can be successively updated by the following function:

$$\widehat{Q}\left(x\right)={\left({\overline{g}}^k\right)}^{T}x$$E3

where ${\overline{g}}^k$ denotes a smoothed estimate of the gradient of the expected recourse function at $x$ for iteration $k$. This method can be proven convergent by projection [9] or recursive linearization [10], although the drawback of this method is that it is time-consuming.

The third group mainly consists of primal and dual decomposition methods. The use of primal and dual decomposition methods dates back to Benders decomposition [11]. Van Slyke and Wets [12] first adopted the L-shaped algorithm into the application of Benders decomposition to two-stage stochastic programs. Pereira and Pinto [13] proposed the stochastic dual dynamic programming (SDDP) method, which has been widely applied in many areas. SDDP uses Benders cuts to compute an outer approximation of a (convex) recourse function, and constructs feasible dynamic programming policies. SDDP has led to numerous related approximation methods those are based on the same logic but seek to improve the approximation procedures by exploiting the underlying structure of the particular applications. These methods consist of use of inexact cuts [14], risk-averse variants [15], embedding SDDP in the scenario tree framework [16]. The convergence of SDDP and related methods has been proven by [17], for linear programs by Girardeau et al. [18].

The fourth group includes separable approximation methods [19, 20]. This type of methods usually replaces the expected recourse function in Eq. (1) with separable approximation functions as follows:

$$\widehat{Q}\left(x\right)=\sum _{i=1}^I{\widehat{Q}}_i\left(x\right)$$E4

If the separable functions ${\widehat{Q}}_i\left(x\right)$ are piecewise linear or linear, we can replace the expected recourse function in Eq. (1) with $\widehat{Q}\left(x\right)$. Then we can solve the problem as a pure network flow problem for network recourse problems, which is polynomial solvable. Thus, it is very efficient. For example, Godfrey and Powell [21] proposed an adaptive piecewise concave approximation (CAVE) algorithm, and the experimental performance of the algorithm shows exceptionally good. However, there was none provable convergent results in their study. In order to provide convergent solutions, Cheung and Powell [19] proposed an approximation algorithm (SHAPE), which uses sequences of strongly convex approximation functions. However, the strongly convex functions require to construct a nonlinear term, and the strongly convex term might damage the pure network structure and need additional computational effort. This chapter intends to introduce an accurate and efficient approximations with the convergence property.

1.3 Contributions of the algorithms

In this chapter, we aim to develop a convergent method that can efficiently approximate the expected recourse function for two-stage stochastic programs. The main contributions are listed as following:

1. We propose a new convergent hybrid learning algorithm to approximate the expected recourse function for two-stage stochastic programs.

2. Through rigorous mathematical analysis, we prove the convergence of the proposed algorithm for general two-stage stochastic programs. The computational results and mathematical analysis both reveals that the algorithm can drop the projection step without sacrificing the convergence for two-stage stochastic programs with network recourse if we can properly construct the initial piecewise linear approximation functions. That means that a pure piecewise linear approximation can be indeed convergent, which is highly consistent with industry practices. This interesting finding answers the open question which has puzzled scholars for more than a decade: why does the piecewise linear approximation work well for two-stage stochastic programs in industry? Our mathematical analysis can provide the first theoretical support.

3. A series of performance analysis has been conducted. The computational results reveal the efficiency of the proposed algorithms and the proposed algorithms are distribution-free. Furthermore, the convergence rate can be affected by the granularity of the initial function ($\delta$). Small granularity usually leads to a high convergence rate. Finally, the computational results also show that the proposed algorithm is very competitive for high dimensional problems.

4. Compared with MAT, the proposed algorithm can collect the information based on knowledge gradients and use it to update the recourse function by learning steps. It can overcome the “curse of dimensionality”. Moreover, it can transfer the problem into a polynomial solvable problem.

The remainder of this chapter is organized as follows. Section 2 presents the description and convergence analysis of the algorithm for general two-stage stochastic programs. The algorithm (without projection steps) for two-stage stochastic programs with network recourse are shown in Section 3. Section 4 demonstrates computational experiments based on an application of the empty container repositioning problem. Section 5 presents the conclusions and outline directions for future research.

2. Description and convergence analysis of ProSHLA for general two-stage stochastic programs

In this section, ProSHLA is first introduced. Subsequently, we analyze the convergence of ProSHLA for general two-stage stochastic programs.

2.1 Description of ProSHLA

To present ProSHLA mathematically, we let, at each iteration $k$,

${\alpha }_k=$ (possibly random) positive step size;

$\overline{Q}\left(x\right)=$ expected recourse function, that is, $E_{\omega }\left[Q\left(x,\omega \right)\right]$;

${\widehat{Q}}^k\left(x\right)=$ a convex differentiable approximation of $\overline{Q}\left(x\right)$;

${\widehat{q}}^k\left(x\right)=$ a subgradient of ${\widehat{Q}}^k\left(x\right)$ at x, that is ${\widehat{q}}^k\left(x\right)$ ∈ ∂${\widehat{Q}}^k\left(x\right)$;

${\overline{g}}^k=$ a smoothed estimate of the gradient of $\overline{Q}\left(x\right)$ at iteration k;

$g^k=$ a stochastic subgradient of $\overline{Q}\left(x\right)$ at $x^k$, that is, g^k ∈ ∂Q(x^k,ω^k + 1);

$H_k=\{{\omega }_1$,${\omega }_2$,…${\omega }_N\}=\text{the history}\mathrm{up}\text{to}\left(\text{and including}\right)\text{iteration}k.$

For a general non-smooth convex function $\widehat{Q}\left(x\right)$, its sub-differential can be defined as follows:

$$\partial \widehat{Q}\left(x\right)=\left\{\widehat{q}\left(x\right)\in {\mathfrak{R}}^n:\widehat{Q}\left(y\right)-\widehat{Q}\left(x\right)\ge \widehat{q}{\left(x\right)}^{\mathrm{T}}\left(y-x\right)\right\}.$$

We combine Eqs. (1), (3), and (4) to form an approximation at iteration k as follows:

$$\mathit{min}{c}_0^{T}x+{\widehat{Q}}^0\left(x\right)+{\left({\overline{g}}^k\right)}^{T}x$$E5

In this study, we approximate the expected recourse function at iteration $k$ via a convex, differentiable approximation ${\widehat{Q}}^0\left(x\right)$ with a linear correction term${\left({\overline{g}}^k\right)}^{T}x$. At each iteration, the linear correction term${\left({\overline{g}}^k\right)}^{T}x$ are introduced to improve the initial approximation ${\widehat{Q}}^0\left(x\right)$. Note that here we use a convex initial approximation function ${\widehat{Q}}^0\left(x\right)$, whereas SHAPE uses strongly convex approximation functions. SHAPE will introduce a nonlinear term in the approximation function to maintain the strong convexity property, and it might destroy the pure network flow problem structure and demands additional computational effort. Moreover, we do not calculate ${\overline{g}}^k$ in the usual manner to obtain stochastic sub-gradients in this study. We use the following form in our model instead:

$$\mathit{min}{c}_0^{T}x+{\widehat{Q}}^k\left(x\right)+{\alpha }_k{\left(g^k-{\widehat{q}}^k\left(x\right)\right)}^{T}x,$$E6

where ${\widehat{Q}}^k\left(x\right)$ is updated as follows:

$${\widehat{Q}}^{k+1}\left(x\right)={\widehat{Q}}^k\left(x\right)+{\alpha }_k{\left(g^k-{\widehat{q}}^k\left(x\right)\right)}^{T}x$$E7

The greatest merit of updating ${\widehat{Q}}^{k+1}\left(x\right)$ in the above way is that it can retain the stochastic sub-gradients $({\widehat{q}}^k\left(x^k\right),{\widehat{q}}^{k-1}$ (x^k − 1),…,${\widehat{q}}^0\left(x^0\right))$ used in the previous iterations. Thus, in iteration $k$, the objective function involves a weighted average of stochastic sub-gradients in the past ($k-1$) iterations. As shown later in Lemma 2, ${\overline{g}}^k$ in Eq. (5) is a linear combination of $g^1,g^2,$…$g^{k-1}$.

Let $P_X:R^n\to X$ be the orthogonal projection onto $X$ [9]. Then, we can obtain a sequence of solutions {$x^k$} using the following procedure (Figure 1).

Figure 1.

The procedure of ProSHLA.

Generally, ProSHLA consists of two-level loops. In the first-level loop, there exists a series of passes, and in the second-level loop, the exists a series of projection steps, which include the step 5 and 6. We first construct an initial bounded and piecewise linear convex approximation function ${\widehat{Q}}^0$ (x) at the beginning of the first pass, then the initial solution x⁰ can be obtained by solving problem (1). A realization of the random quantity ω ∈ Ω can be drawn, and then we can obtain a stochastic sub-gradient of ${\overline{Q}}^0\left(x\right)$ by solving the resulting deterministic problem. Compared with the slope of ${\widehat{Q}}^0$ (x) and the stochastic sub-gradient at x = x⁰, the difference of these two slopes can be used as a linear term to update ${\widehat{Q}}^0$ (x). Subsequently, we can obtain a new solution x^k + 1 using the updated approximation function. If the sub-gradient vectors ${\widehat{q}}^{\overline{m}}\left(x^{k+1}\right)$ of the newly obtained solution x^k + 1 are equal to sub-gradient of solution $x^{\overline{m}}$, that is ${\widehat{q}}^{\overline{m}}\left(x^{\overline{m}}\right)$, the piecewise linear approximation might have jumped into a local optimum. Subsequently, ProSHLA need to jump out from local optimum by implementing projection steps in the second-level loop. If we obtain a new solution x^k + 1 in the second-level loop and the sub-gradient ${\widehat{q}}^{\overline{m}}\left(x^{k+1}\right)$ is different from sub-gradient ${\widehat{q}}^{\overline{m}}\left(x^{\overline{m}}\right)$, then ProSHLA can jump out the second-level loop, and comes to the end of second pass. Thus then, we can repeat the entire process. Finally, ProSHLA will be terminated when the total absolute change in ${\widehat{Q}}^l\left(x\right)$ over a certain number of iterations is low (e.g. ${\sum }_{l=k-M+1}^k\Vert {\widehat{Q}}^l\left(x\right)-{\widehat{Q}}^{l-1}\left(x\right)\mid <\delta$).

Here we point out the main difference between SHAPE and ProSHLA. The most remarkable difference is that ProSHLA uses convex approximation functions while SHAPE uses strongly convex approximation functions. The strongly convexity always maintains a nonlinear term in the approximation function. And this term might destroy the pure network flow structure and causes additional computational effort. To overcome the drawback of the SHAPE, we introduce the projection step in the second-level loop and construct approximation functions. Particularly, the approximation functions in the ProSHLA is NOT strictly convex, while it needs to be strictly convex in SHAPE. Without the projection step in the second-level loop, the ProSHLA might stuck in the corner local –optimum for stochastic linear programs. Thus, ProSHLA can work well for most practical stochastic linear programs, because most of practical stochastic programs are piecewise convex problems.

2.2 Convergence analysis of ProSHLA

Firstly, we demonstrate the convergence theorem of ProSHLA in this subsection. Then, several properties of approximation are presented. Finally, we use these properties to prove the convergence of ProSHLA.

Without loss of generality, the following assumptions are listed.

(A.1)$X\subset {\mathfrak{R}}^n$ is compact and convex.

(A.2)$E_{\omega }Q\left(x_1,\omega \right)$ is convex, finite and continuous on$X$.

(A.3)$g^k$ is bounded such that $‖g^k‖\le c_1$ for each $\omega \in \mathrm{\Omega }$;${\widehat{q}}^k$ is bounded such that $‖{\widehat{q}}^k‖\le c_2$ for each$\omega \in \mathrm{\Omega }$.

(A.4) Piecewise linear function ${\widehat{Q}}^k\left(x\right)$ are convex, implying that

$${\widehat{Q}}^k\left(x_1\right)-{\widehat{Q}}^k\left(y_1\right)\le {\widehat{q}}^k{\left(x_1\right)}^T\left(x_1-y_1\right).$$

(A.5) The stepsize ${\alpha }_k$ are ${\mathcal{H}}_k$ measurable and satisfy

$$0<{\alpha }_k<1,\sum _{k=0}^{\infty }E\left\{{\alpha }_k^2\right\}\le \infty$$

Except for the assumption from (A.1) to (A.5), we also introduce the following assumption to characterize the piecewise linear convex approximation functions.

(A.6) There exists a positive $b$ and a constant $\delta$, such that for any two points $x_1,y_1\in X,$ if$\left|x_1-y_1\right|>\delta$, then $\left|\widehat{q}\left(x_1\right)-\widehat{q}\left(y_1\right)\right|\ge b\left|x_1-y_1\right|$. If there exists $\widehat{q}\left(x_1\right)$ and $\widehat{q}\left(y_1\right)$ such that ${\widehat{q}}^k\left(x_1\right)-{\widehat{q}}^k\left(y_1\right)=0,$ then $\left|x_1-y_1\right|\le \delta$. If $\delta \to 0$, then the function corresponds to a strongly convex function, if $\delta \to \infty$, then the function becomes purely linear.

Given assumption (A.1)–(A.6), we obtain the following theorem of ProSHLA.

Theorem 1. If assumptions (A.1)–(A.6) are satisfied, then the sequence of $\left\{x_1^k\right\}$ generated by algorithm ProSHLA converges almost surely to the optimal solution $x_1^{\ast }\in X^{\ast }$ of problem (1).

In order to prove the Theorem 1, we need to use the following Martingale convergence theorem and three lemmas.

Martingale Convergence Theorem. A sequence of random variables $\left\{W^k\right\}$, which are ${\mathcal{H}}_k$ measurable, is said to be a super-martingale if there exists the sequence of conditional expectations $E\{\left\{W^{k+1}|{\mathcal{H}}_k\right\}$ and satisfies $E\{\left\{W^{k+1}|{\mathcal{H}}_k\right\}\le W^k$.

Theorem 2. (From reference [22]) Let $W^k$ be a positive super-martingale. Then, $W^k$ converges to a finite random variables a.s.

From the above theorem, we can conclude that $W^k$ is a stochastic decreasing analogue essentially.

Based on the convexity property, the optimal solution for problem (8) at iteration $\overline{m}$ can be characterized by the following inequality:

$${\left({\widehat{q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)\right)}^T\left(x_1-x_1^{\overline{m}}\right)\ge 0,\forall x_1\in X$$E8

To obtain Theorem 1, the following three lemmas are required. The first lemma shows that the difference between the solutions of two consecutive update processes will be bounded by the step-size and the stochastic gradient. The second lemma indicates that the approximation ${\widehat{Q}}^k\left(x_1\right)$ is finite. The third lemma shows that $T^k$ (which will be denoted in Lemma 3) is bounded.

Lemma 1. For any two iterations $j\in \left[\overline{m+1},\overline{m+2}\right)$, $i\in \left[\overline{m},\overline{m+1}\right)$, solutions $x_1^j$ and $x_1^i$ obtained by ProSHLA can be characterized by the following inequality:

$${\alpha }_{\overline{m}}{g}^{\overline{m}}\left(x_1^i-x_1^j\right)\le {\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b$$E9

Proof. Consider a special case, where $i$ and $j$ corresponds to two consecutive iterations. Let $i=\overline{m+1}-1$ and $j=\overline{m+1}$. Based on (10), we can obtain that,

$${\left({\widehat{q}}^{\overline{m+1}}\left(x_1^{\overline{m+1}}\right)\right)}^T\left(x_1-x_1^{\overline{m+1}}\right)\ge 0,\forall x_1\in X$$E10

According to the approximation function’s updating rule, we can conclude that,

$${\left({\widehat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}}\right)+{\alpha }_{\overline{m+1}-1}\left(g^{\overline{m+1}-1}-{\widehat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)\right)\right)}^T\left(x_1-x_1^{\overline{m+1}}\right)\ge 0,\forall x_1\in X$$E11

Substituting $x_1$ with $x_1^{\overline{m+1}-1}$ in Eq. (13), we can obtain

$$\begin{array}{l}{\alpha }_{\overline{m+1}-1}{\left(g^{\overline{m+1}-1}-{\hat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)\right)}^T\left(x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}\right)\\ \ge {\hat{q}}^{\overline{m+1}-1}{\left(x_1^{\overline{m+1}}\right)}^T\left(x_1^{\overline{m+1}}-x_1^{\overline{m+1}-1}\right)\end{array}$$E12

If we arrange the above terms, then we can obtain the inequality below:

$$\begin{array}{l}{\alpha }_{\overline{m+1}-1}{\left(g^{\overline{m+1}-1}\right)}^T\left(x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}\right)\\ \ge {\hat{q}}^{\overline{m+1}-1}{\left(x_1^{\overline{m+1}}\right)}^T\left(x_1^{\overline{m+1}}-x_1^{\overline{m+1}-1}\right)-{\alpha }_{\overline{m+1}-1}{\left({\hat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)\right)}^T\left(x_1^{\overline{m+1}}-x_1^{\overline{m+1}-1}\right)\\ =\begin{array}{l}{\left({\hat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}}\right)-{\hat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)\right)}^T\left(x_1^{\overline{m+1}}-x_1^{\overline{m+1}-1}\right)\\ +\left(1-{\alpha }_{\overline{m+1}-1}\right){\left({\hat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)\right)}^T\left(x_1^{\overline{m+1}}-x_1^{\overline{m+1}-1}\right)\end{array}\end{array}$$E13

When iteration $\overline{m+1}-1$ and $\overline{m+1}$ are not in the same update process, then it means that ${\widehat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)\ne {\widehat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}}\right).$ According to assumption (A.6), we can conclude that $\left|{\widehat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)-{\widehat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}}\right)\right|\ge b\left|x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}\right|$.

According to Eqs. (10) and (13), and$0<{\alpha }_{\overline{m+1}-1}<1$, we can obtain

$$\begin{array}{l}{\alpha }_{\overline{m+1}-1}{\left(g^{\overline{m+1}-1}\right)}^T\left(x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}\right)\\ \ge b{‖x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}‖}^2+\left(1-{\alpha }_{\overline{m+1}-1}\right){\left({\hat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)\right)}^T\left(x_1^{\overline{m+1}}-x_1^{\overline{m+1}-1}\right)\\ \ge b{‖x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}‖}^2.\end{array}$$

Applying Schwartz’ inequality, we can get the inequality below:

$$\begin{array}{l}{\alpha }_{\overline{m+1}-1}‖g^{\overline{m+1}-1}‖\bullet ‖x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}‖\ge {\alpha }_{\overline{m+1}-1}{\left(g^{\overline{m+1}-1}\right)}^T\left(x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}\right)\\ \ge b{‖x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}‖}^2\end{array}$$

Therefore, $‖x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}‖\le {\alpha }_{\overline{m+1}-1}\bullet c_1/b$. We can obtain the inequality below:

$${\alpha }_{\overline{m+1}-1}{\left(g^{\overline{m+1}-1}\right)}^T\left(x_1^{\overline{m+1}-1}-x_1^{\overline{m+1}}\right)\le {\left({\alpha }_{\overline{m+1}-1}{c}_1\right)}^2/b$$E14

For any $i\in \left[\overline{m},\overline{m+1}-1]\right.$, $g^i=g^{\overline{m}}=g^{\overline{m+1}-1}$ and${\widehat{q}}^i\left(x_1\right)={\widehat{q}}^{\overline{m}}\left(x_1\right)={\widehat{q}}^{\overline{m+1}-1}\left(x_1\right)$. Thus,

$${\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^i-x_1^{\overline{m+1}}\right)\le {\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b$$E15

For any $i\in \left[\overline{m},\overline{m+1}-1]\right.$ and $j\in \left[\overline{m+1},\overline{m+2}-1]\right.,g^j=g^{\overline{m+1}}=g^{\overline{m+2}-1}$ for any$x_1\in X$. Thus,

$${\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^i-x_1^j\right)\le {\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b$$E16

□

Lemma 2. In iteration $k,$ the approximation function ${\widehat{Q}}^k\left(x_1\right)$ can be written as ${\widehat{Q}}^k\left(x_1\right)={\widehat{Q}}^0\left(x_1\right)+{\left({\overline{g}}^k\right)}^T{x}_1$, where ${\overline{g}}^k$ is a finite vector.

Proof. According to Eq. (5) in Proposition 1, we can conclude that ${\overline{g}}^{k+1}$ is a linear combination of $g^1,g^2,\dots ,g^k$. Since $g^k$ and ${\widehat{q}}^0\left(x_1\right)$ are finite, there will exists a finite and positive vector $\widehat{d}$ such that

$$\widehat{d}\ge \underset{k}{max}\left|g^k-{\widehat{q}}^0\left(x_1^k\right)\right|$$E17

According to Lemma 2 in [19], we can conclude that ${\overline{g}}^{k+1}\ge \widehat{d}$. □

Let$T^k={\widehat{Q}}^k\left(x_1^{\ast }\right)-{\widehat{Q}}^k\left(x_1^k\right)$, where $x_1^{\ast }$ represents the optimal solution. The following Lemma characterizes the difference between $T^{k+1}$ and$T^k$.

Lemma 3. For any two iterations $i\in \left[\overline{m-1},\overline{m}-1]\right.$ and $j\in \left[\overline{m},\overline{m+1}-1]\right.$, $T^i$ and $T^j$ satisfy

$$T^j-T^i\le {\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^i-x_1^j\right)+{\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^i\right).$$E18

Proof. The special case is first considered. Let $i=\overline{m}$ and $j=\overline{m+1}$. By re-writing $x_1^{\ast }-x^{\overline{m+1}}$ as $x_1^{\ast }-x^{\overline{m}}+x^{\overline{m}}-x^{\overline{m+1}}$, we can obtain the following equation:

$$\begin{array}{rcl}{\hat{Q}}^{\overline{m+1}}\left(x_1\right)& =& {\hat{Q}}^{\overline{m+1}-1}\left(x_1\right)+{\alpha }_{\overline{m+1}-1}{\left(g^{\overline{m+1}-1}-{\hat{q}}^{\overline{m+1}-1}\left(x_1^{\overline{m+1}-1}\right)\right)}^T{x}_1\\ & =& {\hat{Q}}^{\overline{m}}\left(x_1\right)+{\alpha }_{\overline{m}}{\left(g^{\overline{m}}-{\hat{q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)\right)}^T{x}_1\end{array}$$

Then,

$$\begin{array}{rcl}{T}^{\overline{m+1}}-T^{\overline{m}}& =& {\hat{Q}}^{\overline{m}}\left(x_1^{\ast }\right)+{\alpha }_{\overline{m}}{\left(g^{\overline{m}}-{\hat{q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)\right)}^T{x}_1^{\ast }-({\hat{Q}}^{\overline{m}}\left(x_1^{\overline{m+1}}\right)\\ & & +{\alpha }_{\overline{m}}{\left(g^{\overline{m}}-{\hat{q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)\right)}^T{x}_1^{\overline{m+1}})-\left({\hat{Q}}^{\overline{m}}\left(x_1^{\ast }\right)-{\hat{Q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)\right)\\ & =& {\alpha }_{\overline{m}}{\left(g^{\overline{m}}-{\hat{q}}^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^{\overline{m}}+x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)+{\hat{Q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)-{\hat{Q}}^{\overline{m}}\left(x_1^{\overline{m+1}}\right)\\ & =& \underset{{\rm I}}{\underbrace{\left({\hat{Q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)-{\hat{Q}}^{\overline{m}}\left(x_1^{\overline{m+1}}\right)-{\alpha }_{\overline{m}}{\left({\hat{q}}^{\overline{m}}\right)}^T\left(x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)\right)}}\\ & & -\underset{\mathit{{\rm I}{\rm I}}}{\underbrace{{\alpha }_{\overline{m}}{\left({\hat{q}}^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^{\overline{m}}\right)}}+\underset{\mathit{{\rm I}{\rm I}{\rm I}}}{\underbrace{{\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)}}+\underset{\mathit{{\rm I}V}}{\underbrace{{\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^{\overline{m}}\right)}}\end{array}$$

Considering each part individually, given that ${\widehat{q}}^{\overline{m}}\in \partial {\widehat{Q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)$, by convexity of ${\widehat{Q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)$, we can obtain

$${\widehat{Q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)-{\widehat{Q}}^{\overline{m}}\left(x_1^{\overline{m+1}}\right)\le {\left({\widehat{q}}^{\overline{m}}\right)}^T\left(x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)$$E19

Thus, the following expression is applicable.

$$\begin{array}{rcl}{\hat{Q}}^{\overline{m}}\left(x_1^{\overline{m}}\right)-{\hat{Q}}^{\overline{m}}\left(x_1^{\overline{m+1}}\right)& \le & {\left({\hat{q}}^{\overline{m}}\right)}^T\left(x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)\\ & =& \left(1-{\alpha }_{\overline{m}}\right){\left({\hat{q}}^{\overline{m}}\right)}^T\left(x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)+{\alpha }_{\overline{m}}{\left({\hat{q}}^{\overline{m}}\right)}^T\left(x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)\end{array}$$E20

Given Eq. (10) and$0<{\alpha }_{\overline{m}}<1$, we know that $\left({\rm I}\right)\le 0$. Additionally, from Eq. (10) and$0<{\alpha }_{\overline{m}}<1$, we know that $\left(\mathit{{\rm I}{\rm I}}\right)\ge 0$. Thus,$T^{\overline{m+1}}-T^{\overline{m}}\le {\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)+{\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^{\overline{m}}\right)$.

For any$i\in \left[\overline{m},\overline{m+1}-1]\right.$, $g^i=g^{\overline{m}}=g^{\overline{m+1}-1}$ and${\widehat{q}}^i\left(x_1\right)={\widehat{q}}^{\overline{m}}\left(x_1\right)={\widehat{q}}^{\overline{m+1}-1}\left(x_1\right)$ for any$x_1\in X$. Therefore,

$$T^{\overline{m+1}}-T^i\le {\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^i-x_1^{\overline{m+1}}\right)+{\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^i\right)$$E21

For any $i\in \left[\overline{m},\overline{m+1}-1]\right.$ and$j\in \left[\overline{m+1},\overline{m+2}-1]\right.,{\widehat{Q}}^j={\widehat{Q}}^{\overline{m+1}}={\widehat{Q}}^{\overline{m+2}-1}$ for any$x_1\in X$. Therefore,

$$T^j-T^i\le {\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^i-x_1^j\right)+{\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^i\right)$$E22

□

To the proof of Theorem 1, we here consider two scenarios. For the first scenario, ProSHLA does not stop in a given update process. Thus, any update process exhibits finite iterations before the algorithm stops, which means $\overline{m+1}-\overline{m}<M$ for any $m$ ($M$ represents a large number). For the second scenario, ProSHLA might terminate in a given update process. In the following text, Theorem 1 is proven for each scenario.

Scenario 1: ProSHLA does not stop in a given update process.

In the first scenario, a subsequence of $\left\{x_1^k\right\},\left\{x_1^{\overline{m}}\right\}$ are considered. We will prove that the subsequence $\left\{x_1^{\overline{m}}\right\}$ converges to the true optimal $x_1^{\ast }.$ According to the definition of $g^k\in \partial Q\left(x_1^k,{\omega }^{k+1}\right)$, we can obtain the following inequality

$${\left(g^k\right)}^T\left(x_1^{\ast }-x_1^k\right)\le Q\left(x_1^{\ast },{\omega }^{k+1}\right)-Q\left(x_1^k,{\omega }^{k+1}\right)$$E23

where $Q\left(x_1,{\omega }^{k+1}\right)$ represents the operational cost function given the outcome${\omega }^{k+1}$.

According to Lemma 1, we can obtain the following inquality:

$${\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^i-x_1^j\right)\le {\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b$$E24

On the basis of Lemma 3, the difference $T^{\overline{m+1}}-T^{\overline{m}}$ can be described as follows:

$$\begin{array}{rcl}{T}^{\overline{m+1}}-T^{\overline{m}}& \le & {\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\overline{m}}-x_1^{\overline{m+1}}\right)+{\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^{\overline{m}}\right)\\ & \le & -{\alpha }_{\overline{m}}\left(Q\left(x_1^{\overline{m}},{\omega }^{\overline{m+1}}\right)-Q\left(x_1^{\ast },{\omega }^{\overline{m+1}}\right)\right)+{\alpha }_{\overline{m}}{\left(g^{\overline{m}}\right)}^T\left(x_1^{\ast }-x_1^{\overline{m}}\right)\\ & \le & -{\alpha }_{\overline{m}}\left(Q\left(x_1^{\overline{m}},{\omega }^{\overline{m+1}}\right)-Q\left(x_1^{\ast },{\omega }^{\overline{m+1}}\right)\right)+{\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b\end{array}$$E25

Conditional expectation of Eq. (27) with respect to ${\mathcal{H}}_k$ can be taken on both side and then we can obtain

$$E\left(T^{\overline{m+1}}|{\mathcal{H}}_{\overline{m}}\right)\le T^{\overline{m}}-{\alpha }_{\overline{m}}\left(\overline{Q}\left(x_1^{\overline{m}}\right)-\overline{Q}\left(x_1^{\ast }\right)\right)+{\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b$$

where $\overline{Q}\left(x_1\right)$ represents the expected recourse function, that is $E_{\omega }Q\left(x_1,\omega \right).$ Given the conditioning on ${\mathcal{H}}_k,T^{\overline{m}}$,${\alpha }_{\overline{m}}$ and $x_1^{\overline{m}}$ on the right-hand side are deterministic. The conditioning ${\mathcal{H}}_k$ cannot provide any information on${\omega }^{\overline{m+1}}$. Hence, we replace $Q\left(x_1,{\omega }^{\overline{m+1}}\right)$ (for $x_1=x_1^k$ and $x_1=x_1^{\ast }$) with its expectation$\overline{Q}\left(x_1\right)$. Given that ${\alpha }_{\overline{m}}\left(\overline{Q}\left(x_1^{\overline{m}}\right)-\overline{Q}\left(x_1^{\ast }\right)\right)\ge 0$, the sequence

$$W^{\overline{m}}=T^{\overline{m}}+{\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b$$E26

is a positive supermartingale. Theorem 2 implies the almost sure convergence of$W^{\overline{m}}$. Hence,

$$T^{\overline{m}}\to T^{\ast }a.s.$$E27

We perform the summation of Eq. (27) from 0 to $\overline{M}$ and obtain the following inequality:

$$T^{\overline{M}}-T^0\le -\sum _{\overline{m}=0}^{\overline{M}}{\alpha }_{\overline{m}}\left(Q\left(x_1^{\overline{m}},{\omega }^{\overline{m+1}}\right)-Q\left(x_1^{\ast },{\omega }^{\overline{m+1}}\right)\right)+\sum _{\overline{m}=0}^{\overline{M}}{\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b$$E28

We take the expectation of both sides. We take the conditional expectation with respect to ${\mathcal{H}}_{\overline{m}}$ and then over all ${\mathcal{H}}_{\overline{m}}$ for the first term on the right-hand side.

$$\begin{array}{l}E\left(T^{\overline{M+1}}-T^0\right)\\ \le -\sum _{\overline{m}=0}^{\overline{M}}E\left\{E\left\{{\alpha }_{\overline{m}}\left(Q\left(x_1^{\overline{m}},{\omega }^{\overline{m+1}}\right)\right)\right\}\right\}\left.\left.\left.-Q\left(x_1^{\ast },{\omega }^{\overline{m+1}}\right)\right\},{\mathcal{H}}_{\overline{m}}\right\}\right\}+E\left\{\sum _{\overline{m}=0}^{\overline{M}}\frac{{\left({\alpha }_{\overline{m}}{c}_1\right)}^2}{b}\right\}\\ \le -\sum _{\overline{m}=0}^{\overline{M}}E\left\{{\alpha }_{\overline{m}}\left(Q\left(x_1^{\overline{m}},{\omega }^{\overline{m+1}}\right)-Q\left(x_1^{\ast },{\omega }^{\overline{m+1}}\right)\right)|{\mathcal{H}}_{\overline{m}}\right\}\\ +{\left(c_1\right)}^2/b\sum _{\overline{m}=0}^{\overline{M}}E\left\{{\alpha }_{\overline{m}}^2\right\}\end{array}$$

We take the limit as $\overline{M}\to \infty$ and use the finiteness of $T^{\overline{M}}$ and $\sum _{\overline{m}=0}^{\overline{M}}E\left\{{\alpha }_{\overline{m}}^2\right\}$ to obtain

$$\sum _{\overline{m}=0}^{\overline{M}}E\left\{{\alpha }_{\overline{m}}\left(Q\left(x_1^{\overline{m}},{\omega }^{\overline{m+1}}\right)-Q\left(x_1^{\ast },{\omega }^{\overline{m+1}}\right)\right)|{\mathcal{H}}_{\overline{m}}\right\}<\infty$$E29

Given that $Q\left(x_1^{\overline{m}},{\omega }^{\overline{m+1}}\right)-Q\left(x_1^{\ast },{\omega }^{\overline{m+1}}\right)\ge 0$ and$\sum _{\overline{m}=0}^{\overline{M}}E\left\{{\alpha }_{\overline{m}}^2\right\}=\infty \left(a.s.\right)$, there exists a subsequence $\left\{\overline{m}\right\}$ such that

$$\overline{Q}\left(x_1^{\overline{m}}\right)\to \overline{Q}\left(x_1^{\ast }\right)a.s.$$

By continuity of $\overline{Q}$, the sequence converges. Hence,

$$x_1^{\overline{m}}\to x_1^{\ast }a.s.$$

Subsequently, we construct another subsequence $\left\{x_1^{\overline{m}-1}\right\}$. Based on Eq. (27),

$$\begin{array}{l}E\left(T^{\overline{M+2}-1}-T^{\overline{M+1}-1}\right)\le -\sum _{\overline{m}=0}^{\overline{M}}{\alpha }_{\overline{m}}\left(Q\left(x_1^{\overline{m+1}-1},{\omega }^{\overline{m+2}-1}\right)-Q\left(x_1^{\ast },{\omega }^{\overline{m+2}-1}\right)\right)+{\left({\alpha }_{\overline{m}}{c}_1\right)}^2/b\end{array}$$

Like-wise, the following approximation can be proved:

$$x_1^{\overline{m}-1}\to x_1^{\ast }a.s.$$

By analogic condition, a very general subsequence $\left\{x_1^i\right\},i\in \left[\overline{m},\overline{m+1}-1\right]$ will almost surely converge to $x_1^{\ast }$. Here, we term this type of subsequence $X_s$.

In the procedure of ProSHLA, the number of all update iterations is finite. Thus, for any subsequence of $\left\{x_1^k\right\}$, we can obtain a subsequence that always belongs to $X_s$. Then, the following conclusion can be obtained:

$${\mathrm{x}}_1^{\mathrm{k}}\to {\mathrm{x}}_1^{\ast }\mathrm{a}.\mathrm{s}.$$

Scenario 2: ProSHLA halts in a given update process.

For the second scenario, ProSHLA halts in a projection procedure which generates a convergent sequence.

Hence, the conclusion of Theorem 1 can be finally obtained. □.

The above analytical processes demonstrate the convergence property of ProSHLA. According to the above results, we require the function $\widehat{Q}\left(x\right)$ to be piecewise linear convex. However, for practitioners are interested in a practically scenarios, in which they usually use separable functions to approximate the expected recourse function for stochastic programs with network recourse. Based on Eq. (4), if the separable functions are piecewise linear or purely linear, then practitioners can easily solve this network recourse problem, because a pure network flow problem is polynomial solvable. In the following section, we will discuss this special practice scenario.

3. Application for two-stage stochastic programs with network recourse using separable piecewise linear functions

In this section, we will discuss the scenario where $\widehat{Q}\left(x\right)$ is separable for two-stage stochastic programs with network recourse. For this scenario, we can simplify implement ProSHLA without projection step. We denote this simplified version as the Stochastic Hybrid Learning Algorithm (SHLA), which is described as follows (Figure 2).

Figure 2.

Description of SHLA.

Essentially, SHLA is not convergent. However, if it is applied to two-stage stochastic programs with network recourse, SHLA will enjoys several merits as follows: (1) the solution of Q(x,ω) is naturally integer; (2) at each iteration, problem Q(x,ω) is simple network flow problem that can be solved by polynomial algorithm.

Here, if we use separable functions, then assumption (A.6) can be satisfied by the following artificially expression:

$${\widehat{q}}^0\left(x_i\right)<{\widehat{q}}^0\left(x_i+\delta \right)$$

Note that for both ProSHLA and SHLA, it allows to choose initial approximation function with different value of $\delta$ flexibly. Thus, if $\delta$ is set to be 1 for any $i$, then we can guarantee the following expression:

$${\widehat{q}}_i^0\left(x_i\right)<{\widehat{q}}_i^0\left(x_i+\delta \right).$$EA.7

Then, we can reach Theorem 3 below.

Theorem 3. If (A.7) is satisfied, SHLA is always convergent for two-stage stochastic programs problem with network recourse.

Proof. For any $x,y\in X,$if there are unequal, we can obtain the following expression according to (A.7).

$${\widehat{q}}^k\left(x\right)\ne {\widehat{q}}^k\left(y\right)$$

Thus,

$$\mid {\widehat{q}}^k\left(x\right)-{\widehat{q}}^k\left(y\right)\mid >0$$

Hence, if we set $\delta =1$ and apply ProSHLA for two-stage stochastic programs with network recourse, then ProSHLA can drop the projection step because ${\widehat{q}}^k\left(x\right)$ and ${\widehat{q}}^k\left(y\right)$ are always unequal. In this situation, ProSHLA is equivalent to SHLA, so SHLA is convergent.

According to the above analysis, we have provided first theoretical convergence support for SHLA-type algorithms which are widely used in numerous applications as mentioned in introduction part. Compared with SHAPE, SHLA does not contain any nonlinear terms so that it can be very efficient. Besides, SHLA can automatically maintain the convexity of the approximation function if the initial piecewise linear functions are properly constructed.

4. Experimental results for performance analysis

In this section, we use two experimental designs to evaluate the performance of the algorithms: (1) An empty container repositioning problem which arises in the context of two-stage stochastic programs with network recourse; and (2) a high dimensional resource allocation problem as an extension experiment. In this section, the empty container repositioning problem is first introduced and then we present the efficiency of ProSHLA and SHLA. Sub-sequentially, we present the convergence of ProSHLA and SHLA, and examine how δ affects convergence performance, and compare the performance under different distributions of random demands. Finally, an extension experiment on a high dimensional resource allocation problem is conduncted to evaluate the efficiency our algorithms.

4.1 Problem generator for the empty container repositioning problem

In this subsection, we test our algorithms in an empty container repositioning problem faced by a major Chinese freight forwarder, who need to manage their numerous empty container in a port network which is located in Pearl River Delta in a fixed route from [23]. The port network contains several hubs (large ports) and spokes (small ports). And the demand of empty container is usually uncertain. When the forwarder need to decide the quantity of empty container to ship from one port to another, they did not know the exact demand of container in the future [24, 25]. Thus, we can formulate the problem as a two-stage stochastic programs with network recourse. Before we formally introduce the problem, we present the following notations.

$L$	=	set of ports;
$d_{\mathit{ij}}$	=	demand from port $i$ to port$j$ in stage 1;
$D_{\mathit{ij}}$	=	demand from port $i$ to port$j$ in stage 2;
$s_i$	=	initial number of empty containers at port $i$;
$S_i$	=	number of empty containers at port $i$ in stage 2;
$c_{\mathit{ij}}$	=	cost for moving an empty container from port $i$ to port$j$;
$r_{\mathit{ij}}$	=	profit for moving a laden container from port $i$ to port$j$;
$x_{\mathit{ij}}$	=	number of laden containers shipped from port $i$ to port$j$ in stage 1;
$y_{\mathit{ij}}$	=	number of empty containers shipped from port $i$ to port$j$ in stage 1;

Then, the problem can be formulated as follows:

$$\mathit{min}\sum _{i\in L}\sum _{j\in L}\left\{-r_{\mathit{ij}}{x}_{\mathit{ij}}+c_{\mathit{ij}}{y}_{\mathit{ij}}\right\}+E_{\omega }\left[Q\left(x,\omega \right)\right],$$E30

s.t.,

$$\sum _{j\in L}\left\{x_{\mathit{ij}}+y_{\mathit{ij}}\right\}=s_i,\forall i\in L$$E31

$$\sum _{i\in L}\left\{x_{\mathit{ij}}+y_{\mathit{ij}}\right\}=s_j,\forall j\in L$$E32

$$y_{\mathit{ij}}\ge 0,\forall i,j\in L$$E33

where the recourse function $Q\left(x,\omega \right)$ is given as follows:

$$Q\left(x,\omega \right)=\mathit{min}\sum _{i\in L}\sum _{j\in L}\left\{-r_{\mathit{ij}}{x}_{\mathit{ij}}\left(\omega \right)+c_{\mathit{ij}}{y}_{\mathit{ij}}\left(\omega \right)\right\}$$E34

s.t.,

$$\sum _{j\in L}\left\{x_{\mathit{ij}}\left(\omega \right)+y_{\mathit{ij}}\left(\omega \right)\right\}=s_i,\forall i\in L$$E35

$$\sum _{i\in L}\left\{x_{\mathit{ij}}\left(\omega \right)+y_{\mathit{ij}}\left(\omega \right)\right\}=s_j,\forall j\in L$$E36

$$y_{\mathit{ij}}\left(\omega \right)\ge 0,\forall i,j\in L$$E37

In order to evaluate the algorithm, a set of problem instances are created. In this study, the problem generator creates ports in L in a 100-mile by 100-mile rectangle. We simply use the Euclidean distance between each pair of ports as the corresponding travel distance. We set the holding cost for a demand to 15 cents per time instance. We set the net profit for a demand to 500 cents per mile. The empty cost is set to 40 cents. The demand $D_{\mathit{ij}}$ between locations i and j is set as follows:

$$D_{\mathit{ij}}={\mathit{\text{out}}}_j·{\mathit{\text{in}}}_i·v,$$

where

${\mathit{\text{out}}}_j$= outbound potential for port $j$;

${\mathit{\text{in}}}_i$ = inbound potential for port $i$;

$v$= random variable.

The outbound and inbound potentials for each port represent the capability of the location to generate outbound demand or attract inbound containers. In the generator, We draw the inbound potential, ${\mathit{\text{in}}}_i$, for port i between 0.2 and 1.8 uniformly, while the corresponding ${\mathit{\text{out}}}_j$ is set as ${\mathit{\text{out}}}_j=2-{\mathit{\text{in}}}_i$ . The reason for this setting is that in real-world regions, large inbound flows port usually exhibits small outbound flows. We also include a random number $v$ with mean 30, that is, the typical daily demand between each pair of locations to capture the randomness in demand. In order to test the performance of the algorithms under different distributions, we also evaluate the performance under exponential, normal and uniform distribution. We set the stepsize${\alpha }_k$ to 1/k.

We solve a deterministic network flow problem to construct an initial piecewise linear functions as described in [1], and we replace the random demand by their mean values in the deterministic problem. Then, we can obtain $\overline{S}=\left\{{\overline{S}}_1,{\overline{S}}_2,\dots ,{\overline{S}}_n\right\}$. For each$i\in L,$ we generate the initial approximation function ${\widehat{Q}}_i^0\left(x\right)=c{\left(x-{\overline{S}}_i\right)}^2,x=0,\delta ,\dots ,\mathit{k\delta },\dots .\mathit{K\delta },$ where c is a positive parameter and x ∈ [0,$\mathit{K\delta }$]. In the projection step, a least-squares problem is solved as following:

$$x^{k+1}=\text{argmin}{\left(x^{k+1}-\left(x^k+{\alpha }_k{g}^k\left(x^k\right)\right)\right)}^2,x^{k+1}\in X.$$

4.2 Effectiveness and efficiency performance

To test the efficiency of the algorithm, we use a myopic algorithm, a posterior bound (PB), the L-shaped algorithm [12] and the inexact cut algorithm [15] as benchmarks. The myopic algorithm simply solves a static deterministic assignment problem at the current stage while ignoring uncertainties in the second stage. It is necessary to solve a deterministic network flow problem with all realized demands to obtain PB. Note that such a posterior optimization involves no uncertainty since decisions are allowed to anticipate future demand. Thus, the cost of PB is the lowest and normally unreachable. As for the L-shaped algorithm and the inexact cut algorithm, a group of linear programming problems with valid cuts should be solved.

We use 8 instances, in which the number of empty containers is ranged from 400 to 3200, and the corresponding number of ports is ranged from 5 to 40. For each instance, 2000 samples are implemented and we obtain the solutions of the myopic algorithm and the sample means of PB, the inexact cut algorithm, the L-shaped algorithm, SHLA and ProSHLA. For SHLA, two classes of initial functions with $\delta$ = 1 and $\delta$ = 2 are selected, whereas we select $\delta$ = 2 for ProSHLA.

We show the experiment results on total cost in Table 1. In Table 1, column 1 presents the number of ports, and column 2 shows the number of the empty containers. The PB bounds are contained in column 3. Columns 4–9 contain the solutions achieved by the myopic algorithm, the L-shaped algorithm, the inexact cut algorithm, SHLA-1, SHLA-2 and ProSHLA, respectively. From the table, it clearly demonstrates that the inexact cut algorithm, the L-shaped algorithm, SHLA, and ProSHLA can achieve optimal or very-near-optimal solutions, which are closer to the PB (lowest) bounds than those of the myopic method. Moreover, the solutions of the L-shaped algorithm are the best solutions known, which are slightly better than the inexact cut algorithm because the latter produces valid cuts that are inexact in the sense that they are not as constraining as optimality cuts in the Lshaped algorithm. In addition, the performance of SHLA ($\delta$ = 1) outperforms that of SHLA ($\delta$ = 2) and ProSHLA ($\delta$ = 2), the reason is that small $\delta$ can lead to good performance. A specific discussion with impact of $\delta$ will be demonstrated later. The performance of ProSHLA ($\delta$ = 2) is slight better than SHLA ($\delta$ = 2). Because the projection steps in ProSHLA help improve the solution. Considering the speed of convergence is quite important in practical problems, we will focus on the computational time for different algorithms, which is shown in the Table 2 below.

$N_L$	$N_R$	Total cost (dollars)
		PB	Myopic	L-shaped	Inexact cut	SHLA-1	SHLA-2	ProSHLA
5	400	−28,551,302	−27,761,331	−28,551,274	−28,551,227	−28,464,248	−28,463,941	−28,464,002
10	800	−59,397,423	−58,432,159	−59,396,702	−59,396,319	−59,338,653	−59,338,190	−59,338,341
15	1200	−98,451,193	−93,188,576	−98,449,868	−98,449,427	−98,257,484	−98,257,244	−98,257,395
20	1600	−147,390,005	−141,062,187	−147,388,360	−147,387,532	−147,269,239	−147,269,212	−147,269,213
25	2000	−185,223,883	−180,875,614	−185,220,423	−185,219,663	−185,092,289	−185,092,113	−185,092,143
30	2400	−234,005,740	−226,978,090	−234,004,427	−234,003,513	−233,401,849	−233,401,516	−233,401,658
35	2800	−266,375,728	−260,966,356	−266,375,468	−266,374,497	−266,337,862	−266,337,649	−266,337,660
40	3200	−304,910,355	−293,882,881	−304,908,597	−304,906,829	−304,907,153	−304,906,829	−304,906,962

Table 1.

Total cost for SHLA and ProSHLA.

NR	NL	Computation time (s)
		Inexact cut	L-shaped	ProSHLA	SHLA-1	SHLA-2
400	5	76	153	43	28	28
800	10	382	535	148	90	95
1200	15	598	1140	332	224	217
1600	20	1084	2277	628	417	420
2000	25	1843	4190	1107	676	790
2400	30	2338	5491	1559	1112	1066
2800	35	4249	8075	2341	1539	1531
3200	40	8154	18,636	5307	3010	3428

Table 2.

Computational time of ProSHLA and SHLA.

As shown in Table 2, ProSHLA and SHLA are more efficient than the inexact cut algorithm and the L-shaped algorithm because ProSHLA and SHLA can utilize the network structure while using the stochastic sub-gradient to approximate the recourse function. From the table, we find that the inexact cut algorithm and L-shaped algorithm are time-consuming, the reason is that here are 2000 samples, and it corresponds to a very large number of cuts for the inexact cut algorithm and L-shaped algorithm. It can also be observed that the computational time of the inexact cut algorithm is smaller than that of the L-shaped algorithm, and the reason is that the optimality cut in L-shaped is more than the valid cuts in the inexact cut algorithm. Moreover, the computational time of SHLA ($\delta$= 1) is almost equal to that of ($\delta$ = 2), which reveals that the computational time canot be affect by the choice of $\delta .$ In contrary, ProSHLA($\delta$ = 2) requires more computational time than SHLA ($\delta$ = 2), and the reason is that the projection step in ProSHLA are time-consuming. In the following text, we focus on the convergence performance of SHLA and ProSHLA. Thus, only the results of the myopic algorithm, PB, ProSHLA and SHLA are demonstrated, and we use the solutions of the myopic algorithm and PB as the upper and lower bounds, respectively.

4.3 Analysis of convergence performance

In this subsection, a set of experiments are conducted to evaluate the convergence performance of SHLA and ProSHLA, and we choose the second instance (N_R = 800 and N_L = 10) as the experimental illustration. The range of the sample number is set from 20 to 640, and we record the result of each combination of N_R and N_L at each iteration. We can seem from Figure 3 that the convergence rate of SHLA-1 and ProSHLA is remarkably high.

Figure 3.

Convergence rate of ProSHLA and SHLA.

To further evaluate how $\delta$ affects the algorithm’s convergence performance, a set of computational experiments are conducted. We increase $\delta$ from 1 to 16 and the number of samples from 20 to 640. Here are many combination of $\delta$ and the number of samples. We record the sample means of the solutions of SHLA and ProSHLA, PB and the myopic method for each combination. We demonstrate the 3D plots of the solution in Figures 4 and 5. As in Figure 4, the layers of ProSHLA and SHLA are extremely close to the PB layer, and this implies the ProSHLA and SHLA are convergent rapidly for various $\delta .$ Furthermore, it can been seem that the performance of ProSHLA can slightly exceeds that of SHLA. In order to further investigate the difference between SHLA and ProSHLA, we demonstrate the performance of ProSHLA and SHLA separately in Figure 5 (without the myopic algorithm and PB). As described in Figure 5, the choice of $\delta$ can affect the performance of ProSHLA and SHLA, and a small $\delta$ usually leads to a good solution.

Figure 4.

Gaps to PB for various $\delta .$

Figure 5.

Comparison of ProSHLA and SHLA for various $\delta$.

We provides more details on the convergence performance of ProSHLA and SHLA for various $\delta$ in Table 3 below, which clearly demonstrates that in conjunction with the small $\delta$, the performance of SHLA and ProSHLA is close to that of PB.

Total cost—% gap to PB (computational time in seconds)
$\delta$	PB	Myopic	SHLA	ProSHLA
1	−59,397,423	1.6251%	0.0989%(909.5)	0.0989% (1441.5)
2	−59,397,423	1.6251%	0.0997%(907.9)	0.0995% (1506.5)
4	−59,397,423	1.6251%	0.1001%(901.8)	0.0997% (1503.6)
6	−59,397,423	1.6251%	0.1005%(911.5)	0.1004% (1553.5)
8	−59,397,423	1.6251%	0.1028%(901.7)	0.1024% (1535.6)
16	−59,397,423	1.6251%	0.1189%(920.3)	0.1150% (1565.4)

Table 3.

Performance under various $\delta$ (no. of samples is 2000).

4.4 An extension experiment on a high dimensional resource allocation problem

Due to the limitation of the container setting, an extension experiment on a higher dimensional problem is considered in this subsection. In this problem, there exists several retailers $R$ and many production facilities (with warehouse) $L$. In stage 1, an amount $x_{\mathit{ij}}$ is moved to a warehouse or retailer or location $j$ from production facility $i$ before the retail demand is realized. When we know the consumer’s demand, then $y_{\mathit{ij}}$ products are moved to retailer location $j$ from production facility $i$. Besides, the type of the consumer’s demand at each location $j$ is different, we denote the type as $t\in {\rm T}$, we set the consumer’s demand of type $t$ at location $j$ as $D_j^t,$ and provide $p_i^t$ unit of type $t$ at production location $i.$ We denote the production capacity of location $i$ by${\mathit{cap}}_i.$ This problem is a non-separable problem.

Subsequently, we formulate the problem as follows:

$$\mathit{min}\sum _{i\in L}\sum _{j\in L\cup R}{c}_{\mathit{ij}}^1{y}_{\mathit{ij}}+E_{\omega }\left[Q\left(x,\omega \right)\right]$$E38

subject to

$$\sum _{j\in L\cup R}{x}_{\mathit{ij}}\le {\mathit{cap}}_i,\forall i\in L$$E39

$$\sum _{i\in L}{x}_{\mathit{ij}}=s_j,\forall j\in L\cup R$$E40

$$x_{\mathit{ij}},s_j\ge 0,\forall i\in L,\forall j\in L\cup R$$E41

where the recourse function$Q\left(x,\omega \right)$ is given as follows:

$$Q\left(x,\omega \right)=\mathit{min}\sum _{i\in L\cup R}\sum _{j\in R}{c}_{\mathit{ij}}^2{y}_{\mathit{ij}}-\sum _{i\in R}\sum _{t\in T}{r}_i^t{p}_i^t$$E42

subject to

$$\sum _{j\in R}{y}_{\mathit{ij}}=s_i,\forall i\in L\cup R$$E43

$$\sum _{i\in L\cup R}{y}_{\mathit{ij}}=\sum _{t\in T}{p}_j^t,\forall i\in R$$E44

$$p_j^t\le D_j^t\left(\omega \right),\forall t\in L\cup R,\forall j\in R,t\in T$$E45

In the first stage, we set $c_{\mathit{ij}}^1=c_0^1+c_1^1{d}_{\mathit{ij}}$, where $d_{\mathit{ij}}$ is the Euclidean distance between locations i and j, and $c_0^1$ is the production cost for each product and $c_1^1$ is the transportation cost per mile. For the second stage costs, we set

$$c_{\mathit{ij}}^2=\left\{\begin{array}{c}{c}_1^2{d}_{\mathit{ij}}\mathit{\text{if}}i\in L\mathit{\text{or}}i=j\\ c_0^2+c_1^2{d}_{\mathit{ij}}\mathit{\text{if}}i\in R\mathit{\text{and}}i\ne j\end{array}\right.$$

$c_1^2$ is the transportation cost per mile in the second stage, and $c_0^2$ represents the fixed charge for moving each product from one retailer location to another retailer location. For one unit of the demand type t occurring in retailer location$i$, a revenue$r_i^t$ will be obtained. Our problem instances differ in the number of products and $\mid L\cup R\mid$, and it determines the dimensionality of the recourse function.

Similarly, we use the inexact cut algorithm [15] and the L-shaped algorithm [12] as benchmarks, and these two algorithms are Benders decomposition based methods. Considering the convergence rate is quite important practically, in this part, our main focus is on the speed of convergence. In order to evaluate the speed of convergence of different methods, each algorithm is implemented for 40, 160, 640, 1200, and 4000 iterations, and a side by side comparison of the algorithms has been made when the number of iteration increases. For the L-shaped and inexact cut algorithms, the number of iterations refer to the number of cuts used to approximate the expected recourse function. For ProSHLA ($\delta$ = 2), the number of iterations refer to the number of demand samples used.

Table 4 below shows the experiment results. In the experiment, the L-shaped algorithm has been used to help find the optimal solution. In the table, the numbers denote the percent deviation between the optimal value and the objective value.

			Number of iterations
N|L∪R|	NP	Algorithm	40	160	640	1200	4000	Sec./iter.
|L ∪ R| = 6	10	ProSHLA; L-shaped	13.26; 1.17	8.61; 0	2.93; 0	2.92; 0	2.73; 0	0.01; 0.07
		Inexact cut	0.96	0	0	0	0	0.05
|L ∪ R| = 10	200	ProSHLA; L-shaped	10.58; 1.85	3.01; 0	0.61; 0	0.29; 0	0.12; 0	0.07; 0.26
		Inexact cut	1.31	0	0	0	0	0.21
|L ∪ R| = 20	400	ProSHLA; L-shaped	7.22; 10.46	1.22; 1.16	0.42; 0	0.23; 0	0.05; 0	0.30; 1.13
		Inexact cut	6.63	0.98	0	0	0	1.04
|L ∪ R| = 40	800	ProSHLA; L-shaped	6.03; 23.57	0.82; 3.23	0.34; 0.31	0.17; 0.02	0.02; 0	0.83; 9.53
		Inexact cut	17.16	2.24	0.13	0.01	0	8.72
|L ∪ R| = 100^∗	2000	ProSHLA; L-shaped	5.68; 44.84	0.78; 14.51	0.15; 1.38	0.04; 0.5	0 0.03	2.68; 36.53
		Inexact cut	29.56	8.14	0.91	0.25	0.03	30.98

Table 4.

Percent error over optimal solution with different algorithms costs.

Note. Figures represent the deviation from the best objective value known.

*Optimal solution not found.

For all problem instances, we use the L-shaped algorithm to find the optimal solution. The numbers in the table represent the percent deviation between the objective value and the optimal value obtained after a certain number of iterations. The computational time per iteration are also listed in Table 4. The computational results on 5 scale of dimensionality instances.

In Table 4 above, column 1 presents the number of the locations, and column 2 shows the number of the products. Column 3 presnets the method that we used in the experiment. The percent deviation from the optimal value are contained in columns 4 to 8. Column 9 lists the computational time per iteration. According to results in Table 4, ProSHLA is able to obtain high quality solutions very efficient for different problem instance, and it can maintain the consistent performance in problem of different sizes, especially for large problems. This performance characteristic makes ProSHLA promising for large-scale application. In comparison with these two Benders decomposition-based methods, ProSHLA is competitive for high dimensional problems. The reason is that separable approximations usually scale much more easily to very high dimensional problems. Note that in the first problem instance, when the number of location is 6 and the number of resource is 10 (the inventory in a location might be 0, 1, 2), the result of ProSHLA seems to be breakdown, because the problem instance in this subsection is non-separable, which may introduce errors when we use the separable approximations to approximate the expected recourse function. However, it will not happen on large problems. As for large problems, the separable approximations are nearly continuous, rather than being just piecewise continuous.

According to the above computational results, ProSHLA is a promising method for two-stage stochastic programs, but more comprehensive numerical work is needed before using it in a particular problem. Owing to its efficient performance and simplicity, ProSHLA is a very promising candidate for high-dimensional problems. Moreover, we can use it as an initialization routine method for high-dimensional stochastic programming problems, and it can exploit high-quality initial feasible solution.

5. Conclusion

In this study, we propose an efficient machine learning algorithm for two-stage stochastic programs. This machine learning algorithm is termed as projected stochastic hybrid learning algorithm, and consists of stochastic sub-gradient and piecewise linear approximation methods. We use the stochastic sub-gradient and sample information to update the piecewise linear approximation on the objective function. Then we introduce a projection step, which implemented the sub-gradient methods, to jump out from a local optimum, so that we can achieve a global optimum. By the innovative projection step, we show the convergent property of the algorithm for general two-stage stochastic programs. Furthermore, for the network recourse problem, our algorithm can drop the projection steps, but still maintains the convergence property. The computational results reveal the efficiency of the proposed algorithms and the proposed algorithms are distribution-free. Furthermore, the convergence rate can be affected by the granularity of the initial function ($\delta$). Small granularity usually leads to a high convergence rate. Finally, the computational results also show that the proposed algorithm is very competitive for high dimensional problems. Compared with MAT, the proposed algorithm can collect the information based on knowledge gradients and use it to update the recourse function by learning steps. It can overcome the “curse of dimensionality”. Moreover, it can transfer the problem into a polynomial solvable problem.

Acknowledgments

This research is partially supported by the National Nature Science Foundation of China (Project No. 71701221) and the Natural Science Foundation of Guangdong Province, China, under Grant number: 2019A1515011127.