A Study of Distributionally Robust Multistage Stochastic Optimization
Jianqiu Huang Kezhuo Zhou Yongpei Guan

		
	
	
	
	 A Study of Distributionally Robust Multistage Stochastic Optimization
 
    Jianqiu Huang, Kezhuo Zhou, and Yongpei Guan
	
 

	Department of Industrial and Systems Engineering

	University of Florida, Gainesville, FL 32611, USA

	Emails: jianqiuhuang@ufl.edu, zhoukezhuo@ufl.edu, and guan@ise.ufl.edu

	
    
=======================================================================================================================================================================================================================================






	
In this paper, we focus on a data-driven risk-averse multistage stochastic programming (RMSP) model considering distributional robustness. We optimize the RMSP over the worst-case distribution within an ambiguity set of  probability distributions constructed directly from historical data samples. The proposed RMSP is intractable due to the multistage nested minimax structure in its objective function, so we reformulate it into a deterministic equivalent that contains a series of convex combination of expectation and conditional value at risk (CVaR), which can be solved by a customized stochastic dual dynamic programming (SDDP) algorithm in this paper. As the size of collected data samples increases to infinity, we show the consistency of the RMSP with distributional robustness to the traditional multistage stochastic programming. In addition, 
to test the computational performance of our proposed model and algorithm, we conduct numerical experiments for a risk-averse hydrothermal scheduling problem, the results of which demonstrate the effectiveness of our RMSP framework.
	
	

	
	Key words: multistage stochastic optimization; data-driven decision making; distributional robustness; hydrothermal scheduling






§ INTRODUCTION

The multistage stochastic program (MSP) has been widely studied in literature for providing multi-period optimal decisions  under uncertainty,  since MSP can be naturally adopted to model various real-life applications with periodical decisions, hydrothermal scheduling <cit.>, power system operations <cit.>, transportation <cit.>, and supply chain planning <cit.>.
For MSP models, the sequential decisions are made depending on the realization of stochastic parameters, which are assumed to follow some known probability distributions. 
Once a decision has been made for the previous period, an observation of the stochastic parameters for the current period  becomes available and then the corresponding decision   will be made considering  future uncertainty. The objective of MSP is to minimize the total expected costs incurred by the decisions over the planning horizon subject to a series of constraints, nonanticipativity constraints and other modeling constraints. Readers are referred to <cit.>, <cit.>, <cit.>, and <cit.> for more detailed  MSP structures and properties. 

 
MSP models are generally computationally intractable, since each decision making depends on not only the past parameter realization, but also decisions made in each previous stage.
To simplify the model, the stochastic parameters are assumed to be discrete random variables and modeled via scenario trees, leading to a deterministic equivalent model of MSP, where the corresponding expectations in the objective function convert to finite sums. 
Two types of decomposition algorithms are proposed to solve MSP, scenario-based methods where the sample of realizations is fixed, and sampling-based methods where the sample of realizations is obtained iteratively.  Scenario-based methods use a small set of realizations from the complete sample space to approximately solve the original program, like diagonal quadratic approximation <cit.>, Lagrangian decomposition <cit.>, L-shaped methods <cit.>, and scenario aggregation methods <cit.>. On the other hand, sampling-based methods iteratively draw a subset of realizations from the complete sample space, where statistical bounds are utilized to create convergence criteria, stochastic dual dynamic programming (SDDP) – a Monte Carlo sampling-based method  <cit.>, 
the stochastic decomposition with extension to the multistage case <cit.>,  and progressive hedging <cit.>. These two types of methods are integrated together in <cit.>.
Moreover, multistage distances are introduced and utilized to solve MSP models approximately <cit.>.


Over the last few decades, risk-averse multistage stochastic programming (RMSP) has been attracting significant attentions, due to its advantages over MSP on modeling certain applications where low probability events have high impact.
The classical MSP focuses on the average behavior, however, there are many real-life applications where low probability events have high impact. 
To handle this issue, various risk measures are designed to remedy the limitation of the traditional expectation operator in the MSP objective function, RMSP uses a convex combination of a expectation and a risk measure <cit.>. Seminal novel works  exploring the risk measure properties for RMSP include coherent risk measures <cit.>, time-consistent risk measures <cit.>, and regularity of risk measures <cit.>. Lagrangian relaxation is first proposed to solve the RMSP models <cit.>, and later come the advanced nested L-shaped decomposition algorithms <cit.>.



Both MSP and RMSP  rely on  known probability distributions, which is generally not practical due to the ambiguity of the probability distributions in real-life applications. Hence, ambiguity sets of probability distributions are proposed to model all possible probability distributions within certain range based on a series of historical data samples. 
Several types of ambiguity sets have been proposed, like moment-based  ambiguity sets <cit.>, divergence-measure-based ambiguity sets <cit.>, and other metric-based ambiguity sets <cit.>. 
Seminal works also construct ambiguity sets based on  the relation between risk measures and robust optimization models <cit.>.



Our contribution in this paper is that  we first propose an RMSP formulation by integrating the general MSP model and the distributional ambiguity sets with  L_∞-norm, which optimizes the total expected costs over the worst-case distribution within the ambiguity set. The ambiguity sets are constructed directly from the historical data samples. We then derive an equivalent reformulation of RMSP,  where the objective function is replaced by a convex combination of an expected cost and a conditional value at risk (CVaR). A significant advantage of our reformulation is that it gets rid of the nested multistage minimax structure from the original objective function, leading to a tractable MSP based on the reference distributions constructed from historical data. We prove that   optimal solutions and objective values of our RMSP with distributional robustness converge to those of risk-netural MSP as the size of data samples increases to infinity. In addition, we customize and analyze the SDDP algorithm to solve our proposed RMSP.  Furthermore, we apply RMSP with distributional robustness to the hydrothermal scheduling problem, and implement corresponding computational experiments whose results verify the convergence of our method.


The following sections of this paper are organized as follows. In Section <ref>, we present our RMSP model with distributional robustness and reformulate it into a tractable equivalent, and show the consistency of our RMSP model. In Section <ref>, we customize the SDDP algorithm to solve our RMSP, followed by a convergence analysis of our SDDP algorithm. 
Next, we apply our RMSP to a risk-averse hydrothermal scheduling problem in Section <ref>, and provide the corresponding computational experiment results. Finally we conclude our research in Section <ref>.










§ RISK-AVERSE MULTISTAGE STOCHASTIC PROGRAM

In this section, we present the reformulation of RMSP based on ambiguity sets with L_∞-norm ^t,  ∀ t = 2,…,T,  and then provide the corresponding convergence analysis. 

At each stage t = 2, 3, ⋯, T, we denote the sample space of the stochastic parameters ξ_t as Ω_t = {ξ_t^1, ξ_t^2, ⋯, ξ_t^R}. Similarly to that of <cit.>, the ambiguity set ^t for possible probability distributions f = (f^1, f^2, ⋯, f^R) can be constructed in a data-driven way with two kinds of representations as follows:

    𝒟^t = { f≥ 0: ||f - f_0 ||_∞≤ d_t,    ∑_r=1^R f^r=1	},

or

    𝒟^t = { f≥ 0: d_1^t ≤ ||f/f_0 ||_∞≤ d_2^t , ∑_r=1^R f^r=1	},

where f_0 represents the reference distribution that can be established from the empirical distribution with historical data, d_t is a tolerance that decreases as the data size N_t increases, and d_1^t/ d_2^t are tolerance parameters that increase/decrease to 1 as the data size N_t grows to infinity. 
The decision rules of d_t, d_1^t, and d_2^t are beyond the scope of this paper and readers are referred to <cit.> for more details in this part. 
It is obvious that representations (<ref>) and (<ref>) can be transformed to each other by selecting proper parameters d_t, d_1^t and d_2^t.
Therefore, we focus on the first approach to construct the distributional ambiguity set 𝒟^t in this paper.
We denote the lower and upper bound of the ambiguity set as f_ℓ = f_0 -d_t and f_u = f_0 +d_t, respectively. 
At each stage t =2,…,T, we let _t^ℓ and _t^u-ℓ represent the probability measures induced by f_ℓ and  f_u -f_ℓ, and we assume P_t^ℓ := ∑_r=1^Rf_ℓ^r ∈ (0,1) and P_t^u := ∑_r=1^Rf_u^r >1  to avoid trivial cases. 
In the following, we reformulate the worst-case expected cost over ^t at each stage t into a convex combination of an expected cost and a CVaR. Finally, we show that both the set of optimal solutions and the objective value of RMSP converge to those of risk-neutral MSP as the data sample size grows to infinity, respectively.




 §.§ Equivalent Reformulation

In this subsection, we develop the equivalent reformulation of the following nested RMSP formulation:  

    min_A_1x_1≥ b_1 c_1x_1 +  sup__2 ∈𝒟^2__2[ min_A_2x_2 ≥ b_2 - B_2x_1 c_2x_2 + ⋯ + sup__T ∈𝒟^T__T[ min_A_Tx_T ≥ b_T -B_Tx_T-1 c_Tx_T  ]  ],

where vectors c_t, b_t, and matrices A_t, B_t are assumed to be stagewise independent random variables forming the stochastic data process (c_t(ξ_t), b_t(ξ_t), A_t(ξ_t), B_t(ξ_t)) for t=2,…, T.


 
	The RMSP has compact feasible set and relatively complete recourse. In addition, the recourse function at each stage is bounded for each decision x. 

Due to the stagewise independence of the data process, formulation (<ref>) can reformulated as a series of dynamic programming equations.
Starting from the last stage T, we define Q_T(x_T-1(ξ_T-1), ξ_T   ) as the optimal value of the last stage program as follows:

    Q_T(x_T-1(ξ_T-1), ξ_T   ) =       min_x_T      c_T(ξ_T)x_T(ξ_T) 
          s.t.       A_T(ξ_T)x_T(ξ_T) ≥ b_T(ξ_T) - B_T(ξ_T)x_T-1(ξ_T-1).



Backward to stage T-1, we have that Q_T-1(x_T-2(ξ_T-2), ξ_T-1   ) is equal to the optimal value of the program

    Q_T-1(x_T-2(ξ_T-2), ξ_T-1) =       min_x_T-1     c_T-1(ξ_T-1)x_T-1(ξ_T-1) +  sup__T∈𝒟^T_ξ_T∼_T[Q_T( x_T-1(ξ_T-1), ξ_T )] 
           s.t.     A_T-1(ξ_T-1)x_T-1(ξ_T-1) ≥ b_T-1(ξ_T-1) - B_T-1(ξ_T-1)x_T-2(ξ_T-2),

and the cost-to-go function _T( x_T-1(ξ_T-1) ) for the worst-case expectation sup__T∈𝒟^T_ξ_T∼_T[  Q_T( x_T-1(ξ_T-1),ξ_T )] can be defined and then reformulated as follows:

    eqn:Qe_T   _T( x_T-1(ξ_T-1) )  ≡sup__T∈𝒟^T_ξ_T∼_T[Q_T( x_T-1(ξ_T-1), ξ_T )] 	
       =  P_T^ℓ_ξ_T∼_T^ℓ[Q_T(x_T-1(ξ_T-1), ξ_T   )] + (1-P_T^ℓ)_(P_T^u-1)/(P_T^u-P_T^ℓ)^ξ_T∼_T^u-ℓ[Q_T( x_T-1(ξ_T-1), ξ_T )]			
        = P_T^ℓ_ξ_T∼_T^ℓ[Q_T(x_T-1(ξ_T-1), ξ_T   )] + inf_ u_T-1∈{  (1-P_T^ℓ)u_T-1(ξ_T-1) 
       + (P_T^u-P_T^ℓ)_ξ_T∼_T^u-ℓ[Q_T(x_T-1(ξ_T-1), ξ_T) - u_T-1(ξ_T-1)]^+    },

where equality (<ref>) holds due to Theorem 6 in <cit.> and equality (<ref>) holds because of the CVaR definition.

Thus, Q_T-1(x_T-2(ξ_T-2), ξ_T-1   ) can be reformulated by substituting Equation (<ref>) into (<ref>).

    min_x_T-1, u_T-1      c_T-1(ξ_T-1)x_T-1(ξ_T-1) + (1-P_T^ℓ )u_T-1(ξ_T-1) + P_T^ℓ_ξ_T∼_T^ℓ[Q_T(x_T-1(ξ_T-1), ξ_T   )] 
       + (P_T^u-P_T^ℓ )_ξ_T∼_T^u-ℓ[Q_T(x_T-1(ξ_T-1), ξ_T) - u_T-1(ξ_T-1)]^+ 
    
    s.t.         A_T-1(ξ_T-1)x_T-1(ξ_T-1) ≥ b_T-1(ξ_T-1) - B_T-1(ξ_T-1)x_T-2(ξ_T-2).


By repeating this process backward, we can write dynamic programming equations for each stage t=2,…,T as

    Q_t(x_t-1(ξ_t-1), ξ_t) =   inf_x_t , u_t{  c_t (ξ_t )x_t (ξ_t ) + (1-P_t+1^ℓ )u_t(ξ_t) + _t+1( x_t (ξ_t ), u_t(ξ_t ) ): 
        A_t (ξ_t )x_t (ξ_t ) ≥ b_t (ξ_t ) - B_t (ξ_t )x_t-1(ξ_t-1) },

where the cost-to-go function _t+1( x_t (ξ_t ), u_t(ξ_t ) ) is defined as   

    _t+1( x_t (ξ_t ), u_t(ξ_t ) ) =    P_t+1^ℓ_ξ_t+1∼_t+1^ℓ[Q_t+1(x_t (ξ_t ), ξ_t+1   )] 
       + (P_t+1^u-P_t+1^ℓ )_ξ_t+1∼_t+1^u-ℓ[Q_t+1(x_t (ξ_t ), ξ_t+1) - u_t(ξ_t) ]^+ ,

with _T+1(·) = 0 and P_T+1^u = P_T+1^ℓ  = 1. 


Finally, the reformulated program at the first stage is described as follows:

    min_x_1, u_1        c_1x_1  +  (1-P_2^ℓ )u_1  + _2( x_1 , u_1 ) 
    
    s.t.         A_1 x_1≥ b_1.



 


 §.§ Convergence Analysis

In this section, we show the consistency of RMSP with distributional robustness by analyzing the convergence property of RMSP as the size of historical data samples increases to infinity. We find that when the size of historical data samples for constructing ambiguity sets at each stage goes to infinity, both the optimal objective value and the set of optimal solutions for RMSP converge to the counterparts of MSP under true while unknown distribution.
We define the following notations for our proof. We let z(0) denote the optimal objective value and U(0) denote the set of optimal solutions for the MSP under  the true distribution. Similarly, we let ẑ(0) denote the optimal objective value and Û(0) denote the set of optimal solutions for the MSP under the reference distribution. We extend the notations to RMSP by denoting ẑ(d_2(N_2),…,d_T(N_T)) as the optimal objective value with data size N_t for each corresponding stage and Û(d_2(N_2),…,d_T(N_T)) as the optimal solutions under the same setting. 

For MSP, the recourse functions are overlapped in a nested structure and only in the last stage do we have a closed-form recourse function for each scenario. Moreover, the number of decision variables in stage t grows in an exponential rate of stage t. Thus, in order to better show our analysis of the convergence property, we denote x_t(ξ_t) as the decision in stage t for observation ξ_t, x_t as the set of decisions in stage t for all observations in stage t, and x as the set of decisions for all stages and scenarios. Besides, we let N=min{N_2,N_3,…,N_T} and d=max{d_2,d_3,…,d_T}. With the above notations, we provide our conclusion on the convergence analysis for RMSP.
 
	As the size of data sample N goes to ∞, ẑ(0)→ z(0). Furthermore, Û(0) converges to U(0), i.e., lim_N →∞sup_x∈Û(0)||x-U(0)||=0.


	Along with the proof in this paper, we will use N→∞ and d→ 0 interchangeably. Note that there exist underlying reference distributions along all stages for each ẑ(0), so we may use the notation ẑ__N(0) when we need to emphasize it. We follow the same notation rule for z(0) and throughout the proof in this paper. We also let ĥ(x,d_2,…,d_T) represent the objective value for RMSP corresponding to solution x and tolerances d_2,…,d_T, and similarly h(x) for MSP with true distribution under solution x.
	
	
    ẑ(0) = ĥ(x̂^*,0) ≤ĥ(x^*,0),

	where x̂^* represents an optimal solution to RMSP under corresponding tolerance, 
	and x^* represents an optimal solution to MSP under true distribution. 
	Taking upper limit with respect to the size of historical data on both sides of Inequality (<ref>), we have
	
    lim sup_N→∞ẑ(0) = lim sup_N→∞ĥ(x̂^*,0) ≤lim sup_N→∞ĥ(x^*,0) = lim_N→∞ĥ(x^*,0) = h(x^*)=z(0), limsup

where the second equality holds
	because of the following reasons: once the decision variables x are fixed, the objective function is a polynomial of parameters f^r_t; the reference distribution for each stage t converges weakly to the true distribution, so the limit of the objective value exists and is consistent with that under true distribution.
	
	Next, we show that {ẑ__N(0)} converges. Under the Assumption <ref> it is easy to see that {ẑ__N(0)} belongs to a bounded set. So if {ẑ__N(0)} does not converge, we can find two subsequences that converge to different values, say z_1 and z_2, i.e., ẑ__N_t(0) → z_1,ẑ__N_s(0) → z_2 and z_1 ≠ z_2. Note that the corresponding optimal solutions  {x̂_N_t^*} and {x̂_N_s^*} are bounded due to the assumption of compact feasible region, there exists subsequences of the two series of optimal solutions that converge respectively. For notation brevity, we still denote the two subsequences as  {x̂_N_t^*} and {x̂_N_s^*}, and x̂_N_t^* → x_1, x̂_N_s^*→ x_2, where x_1 and x_2 are two feasible solution. Then we have
	
		       z_1 = lim_t→∞ ẑ__N_t(0) = lim_t→∞ĥ(x̂_N_t^*,0) = lim_t→∞ĥ(x_1,0) = h(x_1) ≥z(0), 

		       z_2 = lim_s→∞ ẑ__N_s(0) = lim_s→∞ĥ(x̂_N_s^*,0) = lim_s→∞ĥ(x_2,0) = h(x_2) ≥z(0). 

	But we have z_1 ≤ z(0) and  z_2 ≤ z(0) due to (<ref>), so z_1 = z(0) =z_2, which is a contradiction. Thus, {ẑ__N(0)} converges and we have
	
    lim_N→∞ẑ(0) = lim sup_N→∞ẑ(0) = z(0).

	
	Finally, we prove the convergence property of Û(0) to U(0) by contradiction. Supposing that sup_x∈Û(0)||x-U(0)|| does not converge to zero as N grows to infinity, there exists a positive number ϵ_0 and a sequence of optimal solutions {x̂_N_k^*} such that ||x̂_N_k^*-U(0)|| > ϵ_0 for all k. Following the same idea as above, we have a subsequence of the optimal solutions {x̂_N_k^*} that converges. For notation brevity, we still denote the subsequence as  {x̂_N_k^*} and  x̂_N_k^* →x̅, where x̅ is a feasible solution. 
	Then we have lim_k→∞ĥ(x̂_N_k^*,0) = lim_k→∞ĥ(x̅,0) = h(x̅).
	Since {x̂_N_k^*} is a sequence of optimal solutions, we conclude that h(x̅) = z(0) and accordingly x̅∈ U(0). However, as x̂_N_k^* →x̅ and ||x̂_N_k^*-U(0)|| > ϵ_0, we have ||x̅-U(0)|| ≥ϵ_0 > 0, which is a contradiction. This completes the proof.



	For all t≥ 2, as the size of historical data samples N_t increases to ∞, the distance tolerance d_t(N_t) → 0, max_ξ_t ∈Ω_t |f^t(ξ_t)-f_0^t(ξ_t)|→ 0, and ẑ(d_2(N_2),…,d_T(N_T))→ z(0). Furthermore, Û(d_2(N_2),…,d_T(N_T)) converges to U(0), i.e., lim_N_2,…,N_T →∞sup_x∈ U(d_2(N_2),…,d_T(N_T))||x-U(0)||=0.


	To prove that lim_N→∞ẑ(d_2(N_2),…,d_T(N_T)) = z(0), it is enough to show that
	
    lim sup_N→∞ẑ(d_2(N_2),…,d_T(N_T)) ≤ z(0), uppgoal

	and
	
    lim inf_N→∞ẑ(d_2(N_2),…,d_T(N_T)) ≥ z(0).  lowgoal

	First, since ĥ(x̂_N,d^*,d_2(N_2),…,d_T(N_T)) ≥ĥ(x̂_N,d^*,0) ≥ĥ(x̂_N,0^*,0), where x̂_N,d^* and x̂_N,0^* represent the optimal solution to corresponding RMSP under tolerance (d_2(N_2),…,d_T(N_T)) and 0 respectively, by taking lower limit, we have
	
		       lim inf_N→∞ẑ(d_2(N_2),…,d_T(N_T)) 

		    =    lim inf_N→∞ĥ(x̂_N,d^*,d_2(N_2),…,d_T(N_T)) lowlim1

		    ≥    lim inf_N→∞ ĥ(x̂_N,0^*,0) lowlim2

		    =     lim_N→∞ ĥ(x̂_N,0^*,0) lowlim3

		    =     lim_N→∞ ẑ(0)lowlim4 

		    =     z(0), lowlim5
	
	where (<ref>) holds by definition, (<ref>) holds because of the inequality we just provide, and (<ref>) to (<ref>) hold because of Proposition <ref>. Thus, inequality (<ref>) is proved.
	
	Second, we have ĥ(x̂_N,d^*,d_2(N_2),…,d_T(N_T)) ≤ĥ(x,d_2(N_2),…,d_T(N_T)) for any feasible solution x. More specifically, it can be written briefly as 
	
    min  c_1x_1 +  sup__2 ∈𝒟^2__2[ min c_2 x_2 + ⋯ + sup__T ∈𝒟^T__T[ min  c_T x_T ]  ] 
    ≤      c_1x_1 +  sup__2 ∈𝒟^2__2[c_2 x_2 + ⋯ + sup__T ∈𝒟^T__T[  c_T x_T ]  ]

	for all feasible x. Taking upper limit on both sides of inequality (<ref>), we can obtain that 
	 
		       lim sup_N→∞ẑ(d_2(N_2),…,d_T(N_T)) 

		    =     lim sup_N→∞ĥ(x̂_N,d^*,d_2(N_2),…,d_T(N_T)) upplim1 

		    ≤    lim sup_N→∞ { c_1x_1 +  sup__2 ∈𝒟^2__2[c_2 x_2 + ⋯+ sup__T ∈𝒟^T__T[  c_T x_T ]  ] } upplim2 

		    =     lim sup_d_T→0 ⋯lim sup_d_2→0{ c_1x_1 +  sup__2 ∈𝒟^2__2[c_2 x_2 + ⋯+ sup__T ∈𝒟^T__T[  c_T x_T ]  ] }  upplim3 

		    =     lim sup_d_T→0 ⋯lim sup_d_3→0 { c_1x_1 + lim sup_d_2→0 sup__2 ∈𝒟^2__2[c_2 x_2 + ⋯+ sup__T ∈𝒟^T__T[  c_T x_T ]  ] } upplim4 

		    =     lim sup_d_T→0 ⋯lim sup_d_3→0 { c_1x_1 + __2[c_2 x_2 + ⋯+ sup__T ∈𝒟^T__T[  c_T x_T ] ] } upplim5 

		    =     ⋯upplim6 

		    =      c_1x_1 +  __2[c_2 x_2 + ⋯+ __T[  c_T x_T ]  ], upplim7
	
	where (<ref>) holds by definition, (<ref>) holds because of inequality (<ref>), (<ref>) holds because N→∞ is equivalent to d_2(N_2),…,d_T(N_T) → 0, (<ref>) holds because here we only consider the limit when d_2 goes to zero, (<ref>) holds because for stage 2, the corresponding empirical distribution converges weakly to the true distribution and the recourse function for stage 2 is actually a linear expression of f_2^r, which leads to both the existence of the limit and the consistency,  and (<ref>) to (<ref>) holds by taking the limit for stage 3,…, T. Since inequality (<ref>) holds for any feasible solution x, we conclude that 
	 
		       lim sup_N→∞ẑ(d_2(N_2),…,d_T(N_T)) 

		    ≤     minc_1x_1 +  __2[ minc_2 x_2 + ⋯+ __T[ minc_T x_T ]  ] 

		    =     z(0).
	
	Thus, inequality (<ref>) is proved and lim_N→∞ẑ(d_2(N_2),…,d_T(N_T)) = z(0).
	
	The proof for convergence property of optimal solutions follows similar method to that in Proposition <ref> and is thus omitted here.







§ SOLUTION APPROACH

In this section, we first utilize a scenario tree  to model the data process ξ_2, …,ξ_T based on the historical data, and then customize the SDDP approach to solve our RMSP.



We collect the historical data and then generate a finite scenario tree based on these data with two reference distributions ^ℓ and ^u-ℓ. At each stage t=2,…, T,  we denote  N_t as the total number of different scenarios, and let p_ti^ℓ and p_ti^u-ℓ be the corresponding reference probabilities for scenario ξ_ti, i=1,…,N_t.
Therefore, we can rewrite the dynamic programming equations in (<ref>) as follows:

    Q̂_ti(x̅_t-1) = inf_x_t , u_t{    c_t (ξ_ti)x_t (ξ_ti) + (1-P_t+1^ℓ )u_t(ξ_ti) + _t+1( x_t (ξ_ti), u_t(ξ_ti) ): 
        A_t(ξ_ti)x_t (ξ_ti) ≥ b_t (ξ_ti) - B_t (ξ_ti)x̅_t-1},

for i=1,…, N_t, where x̅_t-1 is the current solution obtained from stage t-1, and the corresponding cost-to-go function can be written as follows:

    _t+1( x_t (ξ_tj ), u_t(ξ_tj ) ) = P_t+1^ℓ∑_i=1^N_t+1p_t+1i^ℓ  [Q̂_t+1i(x_t (ξ_tj ))] + (P_t+1^u-P_t+1^ℓ )∑_i=1^N_t+1p_t+1i^u-ℓ[Q̂_t+1i(x_t (ξ_tj )) - u_t(ξ_tj) ]^+.


At the first stage, we have

    min_x_1, u_1        c_1x_1  +  (1-P_2^ℓ )u_1  + _2( x_1 , u_1 ) 
    
    s.t.         A_1 x_1≥ b_1.



In the following, we use variables x_t, u_t instead of x_t (ξ_tj ), u_t(ξ_tj ) for notation brevity, and we denote x̅_t , u̅_t as the current solutions for stage t.
As the cost-to-go functions _t+1( x_t , u_t ) are convex,   we can use supporting hyperplanes to  make lower approximations of the cost-to-go functions. We define 𝔔_t+1(x_t , u_t) as the current approximation of the cost-to-go function at stage t+1 for each t= 1,…,T-1. 

To generate a supporting hyperplane for _t+1( x_t , u_t ) at ( x̅_t , u̅_t ) in (<ref>), we first consider the subdifferential of function [Q̂_t+1i(x_t ) - u_t]^+ at point (  x̅_t , u̅_t ), which is 
	
    ∂ [Q̂_t+1i(x_t ) - u_t]^+  = {[                                 [0,0]             if Q̂_t+1i(x̅_t ) < u̅_t,;           ⋃_g∈∂Q̂_t+1i(x̅_t )  [g,-1]            if Q̂_t+1i(x̅_t ) > u̅_t ,; ⋃_g∈∂Q̂_t+1i(x̅_t ) , ∈ [0,1] [ g, -]             if Q̂_t+1i(x̅_t ) = u̅_t. ].




Applying the chain rule of subdifferentials, the subgradient of _t+1( x_t , u_t ) at ( x̅_t , u̅_t ) is 
	
    [ P_t+1^ℓ∑_i=1^N_t+1p_t+1i^ℓ𝐠_t+1i + (P_t+1^u - P_t+1^ℓ  )∑_i∈𝒥_t+1p_t+1i^u-ℓ𝐠_t+1i,   -(P_t+1^u-P_t+1^ℓ )|𝒥_t+1|  ],

	where 𝐠_t+1i∈∂Q̂_t+1i(x̅_t ), i=1,…,N_t+1 and 𝒥_t+1:={j:  Q̂_t+1i(x̅_t ) > u̅_t}.



	The supporting hyperplane of _t+1( x_t , u_t ) at ( x̅_t , u̅_t ) is constructed as follows.
	
    ℎ_t+1(x_t, u_t) :=    _t+1(x̅_t, u̅_t) +  [P_t+1^ℓ∑_i=1^N_t+1p_t+1i^ℓ𝐠_t+1i + (P_t+1^u-P_t+1^ℓ )∑_i=1^N_t+1p_t+1i^u-ℓ𝐠_t+1i](x_t - x̅_t) 
        -(P_t+1^u-P_t+1^ℓ )|𝒥_t+1| (u_t -u̅_t).



Therefore, we can update the cost-to-go functions by adding the new supporting hyperplane (<ref>), 𝔔_t+1(x_t , u_t) := max{𝔔_t+1(x_t , u_t), ℎ_t+1(x_t, u_t)  }. 





 §.§ Algorithm




In this section, we adopt SDDP to solve our RMSP. SDDP was first proposed in <cit.> to solve MSP and later studied in <cit.> for MSP models where the objective function is a convex combination of expectation and CVaR. There are two parts in SDDP, a forward simulation that generates a simulated solution, and a backward recursion that improves the approximation of the cost-to-go functions at each stage.
The detailed algorithm is provided in Algorithm <ref>. We describe the application of SDDP to our RMSP in the following section.

In the backward recursion steps, we solve the optimization problem (<ref>) with the approximation  𝔔_t+1(x_t , u_t) instead of the cost-to-go function _t+1( x_t , u_t ).

    Q̂_ti(x̅_t-1) = inf_x_t , u_t{    c_t (ξ_ti)x_t (ξ_ti) + (1-P_t+1^ℓ )u_t(ξ_ti) + 𝔔_t+1( x_t (ξ_ti), u_t(ξ_ti) ): 
        A_t(ξ_ti)x_t (ξ_ti) ≥ b_t (ξ_ti) - B_t (ξ_ti)x̅_t-1},

for i=1,…, N_t, where the corresponding approximated expected cost-to-go functions can be updated by adding supporting hyperplanes (<ref>).

At the first stage,

    min_x_1, u_1        c_1x_1  +  (1-P_2^ℓ )u_1  + 𝔔_2( x_1 , u_1 ) 
    
    s.t.         A_1 x_1≥ b_1,

whose optimal objective value is used as the lower bound of our problem.


The forward simulation steps are performed by sampling independent scenarios from the historical data and computing the corresponding optimal value, which will  be used later to calculate the upper bound z̅ = c_1x̂_1 + (1-P_2^ℓ )û_1  + 𝔳̂_2, 
where 
    𝔳̂_t =1M_t∑_i=1^M_t{  P_t+1^ℓ (c_tix̂_ti + 𝔳̂_t+1 i)  + (1- P_t+1^ℓ   )û_ti + (P_t+1^u- P_t+1^ℓ  )[c_tix̂_ti + 𝔳̂_t+1 i -û_ti]^+}, t=2,…,T.







 §.§ Convergence Property

In this subsection, we provide the convergence of the proposed algorithm.


	When the basic optimal solutions are employed in the backward steps, the forward step procedure generates the optimal solution for the risk-averse multistage stochastic program w.p.1 after a sufficiently large number of backward and forward steps of the algorithm.

Since the total number of scenarios generated from the historical data is finite and the forward steps generate the sample scenarios independently, the conditions in Proposition 3.1 in <cit.> hold. Thus, our proposition holds. 




§ RISK-AVERSE HYDROTHERMAL SCHEDULING PROBLEM

In this section, we investigate a risk-averse hydrothermal scheduling problem with distributional robustness by applying our proposed reformulation and algorithm in Sections <ref> – <ref>.

The risk-averse hydrothermal scheduling problem is aiming to determine an optimal operational schedule for the hydrothermal system which minimizes the total expected costs including operation cost, fuel cost, and penalty cost for failing to satisfy the electricity load. The hydrothermal system is a combination of hydroelectric generators and other fuel-costing generators like thermal and nuclear generators. 
Hydroelectric generators utilize stored water in the system reservoirs to provide energy, whose generation is determined by the water inflow. 
The inflow is assumed to be the only uncertainty of this problem, which follows some unknown distribution and we only have historical data of the inflow volume at each stage.
The risk-averse hydrothermal scheduling problem can be modeled as a stochastic dynamic program if the inflow volume of the current stage is predicted at the beginning of that stage. 




 §.§ Problem Formulation

To describe the problem, we let T represent the planning horizon, and 
we denote sets of  thermal generators and reservoirs as 𝒢 and ℛ, respectively.
For each thermal generator g ∈𝒢, we label the maximum (minimum) generation as G^g(G^g), and denote the unit generation cost as c_t^g at time t.
For each reservoir r∈ℛ, we denote the maximum (minimum) reservoir storage as R^r (R^r), and denote the inflow to reservoir r at time t as I^r_t. For convenience, we use energy units, MWh, as the units of R^r (R^r) and I^r_t, which can easily be done by multiplying  constant coefficients of the true storage level and inflow amount, respectively.
Furthermore, we use D_t to represent the total electricity load, and  use c^p_t to represent the penalty cost for each unit of unsatisfied load at time t.
 
Next, we define five continuous decision variables using energy units MWh as follows. We define 
x^g_t as the thermal energy generated by generator g∈𝒢  at time t, and x^p_t as the energy amount failed to satisfy at time t. For each reservoir r∈ℛ, we define  
the hydroelectric energy generation at time t as x^r_t,
the water spillage as s_t^r, and  the storage level as v^r_t for reservoir r∈ℛ at time t, where the initial storage v^r_0 is given as a parameter. The generated hydroelectricity energy is a linear function of the water outflow, which can be represented by the outflow volume times a constant. 

Based on the above notation, we provide the formulation as follows.

    hydro-formulationmin_x_1        (∑_g∈𝒢c_1^g x_1^g +  c_1^p x_1^p) +  sup__2 ∈𝒟^2_ξ_2∼_2[ min_ x_2 (c_2^g(ξ_2)x_2^g(ξ_2)  +  c_2^p(ξ_2)x_2^p(ξ_2))  
       + ⋯ + sup__T ∈𝒟^T_ξ_T∼_T[ min_ x_T (c_T^g(ξ_T)x_T^g(ξ_T)  +  c_T^p(ξ_T)x_T^p(ξ_T))   ]  ]   
    
    s.t.        ∑_r∈ℛx_t^r + ∑_g∈𝒢x_t^g + x_t^p = D_t, ∀ t =1, …,T, 
        v_t^r + x_t^r + s_t^r = v_t-1^r + I_t^r, ∀ t =1, …,T, ∀ r ∈ℛ, 
       G^g ≤ x_t^g ≤G^g, ∀ t =1, …,T, ∀ g ∈𝒢, 
       R^r ≤ v_t^r ≤R^r, ∀ t =1, …,T, ∀ r ∈ℛ, 
        x_t^p, x_t^h, s_t ≥ 0, ∀ t =1, …,T,

where the objective function (<ref>) is to minimize long-term expected thermal generation cost and penalty cost under the worst-case distribution within the ambiguity set. Constraints (<ref>) ensure the electricity load balance. Constraints (<ref>) restrict on the reservoir water balance, where the inflow amount equals to the sum of storage level difference, spillage and water outflow for hydroelectricity generation. Constraints (<ref>) and (<ref>) represent capacities of the thermal generation and the reservoir storage, respectively. Constraints (<ref>) are the nonnegative constraints of decision variables.



 §.§ Experiment Settings


In the following, we perform numerical experiments on the risk-averse hydrothermal scheduling problem with distributional robustness by applying our proposed reformulation and algorithm in Sections <ref>-<ref>.  We randomly generate four cases by assuming the inflow model follows four different classes of true probability distributions, lognormal, truncated normal, Weibull, and exponential distributions. For each case at each stage, we first create and collect N independent random samples from the true distribution to estimate the empirical distribution f_0, and then construct a scenario tree with S scenarios at each stage.
Next, the ambiguity set is constructed as 𝒟 = {f≥ 0: f_0^r -d ≤ f^r ≤ f_0^r+d, ∀ r = 1, …, R, ∑_r=1^Rf^r =1   }, where d = max{z_α/2√(f^r_0(1-f_0^r))/√(N), ∀ r = 1,…,R } and z_α/2 is the α-level z-score (here we select α= 5%, and z_α/2 = 1.96).
Finally, we compare our risk-averse optimal solutions with risk-neutral solutions under perfect information. All experiments were coded in C++ and implemented on a computer node with two AMD Opteron 2378 Quad Core Processors at 2.4GHz and 4GB memory. IBM ILOG CPLEX 12.3 is utilized as the linear programming solver. 

For the parameter setting of Algorithm <ref>, we set the optimality gap as 5%, the iteration limit K as 300, the sample number at each stage M_t, ∀ t =2,…,T as 6. 
For uncertainty set parameter settings,  
we create a scenario tree with 12 different scenarios at each stage, S=12. For the hydrothermal scheduling problem  parameter settings, we consider a 52-week planning horizon for a single reservoir and a single thermal generator. We let the reservoir maximum capacity be R =  10^6 and minimum capacity be R = 10^5. The initial storage is set as v_0 = 5.5 × 10^5. The thermal unit generation costs at each stage are randomly generated which vary from $ 45 / to $ 85 /, and the unit penalty cost is set as $ 1000 /.




 §.§ Computational Results









In the following, we compare our risk-averse solutions and their corresponding risk-neutral solutions with perfect information under various data sample sizes and true distribution settings in Figure <ref>. 
As indicated from Figure <ref>, our risk-averse solutions converges to their corresponding risk-neutral solution as the number collected data samples increases, and the convergence process evolves moderately quickly after collecting 2000 data samples, which numerically proves the convergence of our risk-averse multistage stochastic program with distributional ambiguity to the risk-neutral multistage stochastic program.

We present the solution gaps between the risk-averse (RA) objective value and the risk-neutral (PI) objective value, and the algorithm iteration steps under various data sample size and distribution settings  in Table <ref>. In this table, the column “Gap" represents the gap between RA and PI, = (Z_a - Z_n)/Z_n, where Z_a is the risk-averse solution obtained from Algorithm <ref>, and Z_n is the risk-neutral objective value under the corresponding true distribution. The column “Step" represents the iteration steps when our algorithm stops, where 300 means that the corresponding case ceases due to the predefined iteration step limit. 
As indicated from Table <ref>, the solution gap decreases and the number of required iteration step becomes smaller as the number collected data samples increases, which is coincident with the convergence of the ambiguity set. That is, with more data samples collected, our ambiguity set size shrinks and thus it is faster to solve the risk-averse stochastic program over this ambiguity set. 



We provide the algorithm performance under various distribution settings with 9000 data samples in Figure <ref>. Figures (<ref>) - (<ref>) represent the upper bound and lower bound evolving processes under different true distributions lognormal, normal, exponential and Weibull, respectively. The horizontal axis is the number of performed iteration steps and the vertical axis is the total cost. Since the gap between upper bound and lower bound is significant in the first 30 iteration steps in the Figure <ref>, we provide a closer look at the convergent behavior of Algorithm <ref> in Figure <ref>, by using the log value of total cost in the vertical axis. The algorithm converges quickly in 35 steps and terminates with about 75 steps for each distribution setting. 





§ CONCLUSION

In this paper, we present an equivalent reformulation of RMSP, where we use a convex combination of expectation and CVaR to replace the worst-case expectation. The reformulation prevents repeating min-max patterns in the multistage program. As the size of collected data samples goes to infinity, we show that RMSP converges to the risk-neutral MSP, where the optimal objective value and the set of optimal solutions of RMSP converge to those of risk-neutral MSP. We adopt the SDDP algorithm to solve the reformulated RMSP and provide the convergence property for the algorithm. To test the RMSP computation performance, we implement numerical experiments for the risk-averse hydrothermal scheduling problem under different true distributions, which demonstrate the convergence of our RMSP to risk-neutral MSP as the collected data increase to infinity.

=12pt
plain


