Risk Neutral Reformulation Approach to Risk Averse Stochastic Programming

TL;DR

This paper presents a method to reformulate risk-averse multistage stochastic programming problems into risk-neutral ones by changing the probability measure, demonstrated through an application to power system operation planning.

Contribution

It introduces a novel change-of-measure approach that simplifies risk-averse problems into risk-neutral forms, facilitating easier analysis and solution.

Findings

The change-of-measure approach effectively emphasizes extreme scenarios.

Application to power system planning shows improved problem handling.

Reformulation benefits include computational advantages and better scenario representation.

Abstract

The aim of this paper is to show that in some cases risk averse multistage stochastic programming problems can be reformulated in a form of risk neutral setting. This is achieved by a change of the reference probability measure making ``bad" (extreme) scenarios more frequent. As a numerical example we demonstrate advantages of such change-of-measure approach applied to the Brazilian Interconnected Power System operation planning problem.

Figures

Equations

$$\underset{x\in \mathcal{X}}{\min }\rho [F(x,\omega )],$$

$$\rho (Z):=(1-\lambda )\mathbb{E}[Z]+\lambda {AV@R}_{\alpha }(Z),\lambda \in (0,1),$$

$${AV@R}_{\alpha }(Z)=\underset{u\in \mathbb{R}}{\inf }\left\{u+{\alpha }^{-1}\mathbb{E}[Z-u{]}_{+}\right\},\alpha \in (0,1).$$

$$\rho (Z)=\underset{\zeta \in \mathfrak{A}}{\sup }{\mathbb{E}}_{\zeta }[Z],$$

$${\mathbb{E}}_{\zeta }[Z]:={\int }_{\Omega }Z(\omega )\zeta (\omega )dP(\omega ),$$

$$\underset{x\in \mathcal{X}}{\min }\underset{\zeta \in \mathfrak{A}}{\sup }{\mathbb{E}}_{\zeta }[F_x].$$

$$\underset{\zeta \in \mathfrak{A}}{\max }\underset{x\in \mathcal{X}}{\inf }{\mathbb{E}}_{\zeta }[F_x].$$

$${\mathbb{E}}_{\bar{\zeta }}[F_x]\ge {\mathbb{E}}_{\bar{\zeta }}[F_{\bar{x}}]\ge {\mathbb{E}}_{\zeta }[F_{\bar{x}}],\forall (x,\zeta )\in \mathcal{X}\times \mathfrak{A}.$$

$$\underset{x\in \mathcal{X}}{\inf }{\mathbb{E}}_{\bar{\zeta }}[F_x].$$

$$\bar{\zeta }\in \arg \underset{\zeta \in \mathfrak{A}}{\max }{\mathbb{E}}_{\zeta }[Z].$$

$$\arg \underset{\zeta \in \mathfrak{A}}{\max }{\mathbb{E}}_{\zeta }[Z]=\partial \rho (Z)$$

$$\partial\rho(Z)=(1-\lambda)\{{\mbox{\boldmath$1$}}\}+\lambda\partial({\sf AV@R}_{\alpha})(Z)$$

$$\partial({\sf AV@R}_{\alpha})(Z)=\left\{\zeta:{\mathbb{E}}[\zeta]=1,\begin{array}[]{lll}\zeta(\omega)=\alpha^{-1},&\text{if}&Z(\omega)>{\sf V@R}_{\alpha}(Z),\\ \zeta(\omega)=0,&\text{if}&Z(\omega)<{\sf V@R}_{\alpha}(Z),\\ \zeta(\omega)\in[0,\alpha^{-1}],&\text{if}&Z(\omega)={\sf V@R}_{\alpha}(Z),\end{array}\right.$$

$${V@R}_{\alpha }(Z)=\inf \{t:P(Z\le t)\ge 1-\alpha \}.$$

$$\rho (Z)=\frac{(1-\lambda )}{N}\sum _{i=1}^N{Z}_i+\lambda \bar{u}+\frac{\lambda }{\alpha N}\sum _{i=1}^N[Z_i-\bar{u}{]}_{+},$$

$$\rho (Z)=\frac{(1-\lambda )}{N}\sum _{i=1}^N{Z}_{(i)}+\lambda Z_{(\kappa )}+\frac{\lambda }{\alpha N}\sum _{i=\kappa +1}^N(Z_{(i)}-Z_{(\kappa )})=\sum _{i=1}^N{\mathfrak{q}}_i{Z}_{(i)},$$

$${\mathfrak{q}}_{i}:=\left\{\begin{array}[]{lll}(1-\lambda)/N&{\rm if}&i<\kappa,\\ (1-\lambda)/N+\lambda-\lambda(N-\kappa)/(\alpha N)&{\rm if}&i=\kappa,\\ (1-\lambda)/N+\lambda/(\alpha N)&{\rm if}&i>\kappa.\end{array}\right.$$

$$\min_{x_{1}\in{\cal X}_{1}}f_{1}(x_{1})+\rho_{2|\xi_{[1]}}\left[\inf_{x_{2}\in{\cal X}_{2}(x_{1},\xi_{2})}f_{2}(x_{2},\xi_{2})+\cdots+\rho_{T|\xi_{[T-1]}}\,\Big{[}\inf_{x_{T}\in{\cal X}_{T}(x_{T-1},\xi_{T})}f_{T}(x_{T},\xi_{T})\Big{]}\right],$$

$$\rho_{t|\xi_{[t-1]}}[Z]=(1-\lambda){\mathbb{E}}\left[Z\big{|}\xi_{[t-1]}\right]+\lambda{\sf AV@R}_{\alpha}\left[Z\big{|}\xi_{[t-1]}\right],\;\lambda,\alpha\in(0,1).$$

$$\begin{array}{c}{f}_t(x_t,{\xi }_t):=c_t^{\mathsf{T}}{x}_t,{\mathcal{X}}_1:=\left\{x_1:A_1{x}_1=b_1,x_1\ge 0\right\},\\ {\mathcal{X}}_t(x_{t-1},{\xi }_t):=\left\{x_t:B_t{x}_{t-1}+A_t{x}_t=b_t,x_t\ge 0\right\},t=2,\dots ,T.\end{array}$$

$$\underset{\pi \in \Pi }{\min }\bar{\rho }\left[f_1(x_1)+f_2(x_2({\xi }_{[2]}),{\xi }_2)+\cdots +f_T(x_T({\xi }_{[T]}),{\xi }_T)\right],$$

$$\bar{\rho}[Z]=\rho_{2|\xi_{[1]}}\Big{[}\rho_{3|\xi_{[2]}}\big{[}\cdots\rho_{T|\xi_{[T-1]}}[Z]\big{]}\Big{]}.$$

$$\bar{\rho }[Z]=\underset{Q\in \mathfrak{M}}{\sup }{\mathbb{E}}_Q[Z],$$

$$\underset{\pi \in \Pi }{\min }\underset{Q\in \mathfrak{M}}{\max }{\mathbb{E}}_Q\left[f_1(x_1)+f_2(x_2({\xi }_{[2]}),{\xi }_2)+\cdots +f_T(x_T({\xi }_{[T]}),{\xi }_T)\right].$$

$$\underset{Q\in \mathfrak{M}}{\max }\underset{\pi \in \Pi }{\min }{\mathbb{E}}_Q\left[f_1(x_1)+f_2(x_2({\xi }_{[2]}),{\xi }_2)+\cdots +f_T(x_T({\xi }_{[T]}),{\xi }_T)\right].$$

$$\underset{\pi \in \Pi }{\min }{\mathbb{E}}_{\bar{Q}}\left[f_1(x_1)+f_2(x_2({\xi }_{[2]}),{\xi }_2)+\cdots +f_T(x_T({\xi }_{[T]}),{\xi }_T)\right],$$

$$V_{t}\left(x_{t-1},\xi_{t}\right)=\inf_{x_{t}\in{\cal X}_{t}(x_{t-1},\xi_{t})}\big{\{}f_{t}(x_{t},\xi_{t})+{\cal V}_{t+1}\left(x_{t}\right)\big{\}},$$

$${\cal V}_{t+1}(x_{t}):=\rho_{t+1}\big{(}V_{t+1}(x_{t},\xi_{t+1})\big{)},$$

$$\bar{x}_{t}\in\arg\min_{x_{t}\in{\cal X}_{t}(\bar{x}_{t-1},\xi_{t})}\big{\{}f_{t}(x_{t},\xi_{t})+{\cal V}_{t+1}\left(x_{t}\right)\big{\}}.$$

$${\rho }_{t+1}(Z)=\underset{{\zeta }_{t+1}\in {\mathfrak{A}}_{t+1}}{\sup }{\mathbb{E}}_{{\zeta }_{t+1}}[Z],$$

Full text

Keywords: stochastic programming, risk measures, dynamic equations, saddle point, Stochastic Dual Dynamic Programming algorithm, importance sampling, power system planning problem

1 Introduction

There are many practical problems where one has to make decisions sequentially based on data (observations) available at time of the decision. In the stochastic programming approach the underlying data is modeled as a random process with a specified probability distribution. We can refer to the books [9],[15], and references therein, for a thorough discussion of the Multistage Stochastic Programming (MSP). In the risk neutral formulation of MSP problems, expected value of the total cost is supposed to be optimized (minimized). Of course for a particular realization of the data process the corresponding cost can be quite different from its average. This motivates to consider risk averse approaches where one tries to control high costs by imposing some type of penalty on high cost realizations of the data process.

An axiomatic approach to risk was suggested in the pioneering paper by Artzner et al [1], where the concept of coherent risk measures was introduced. By the Fenchel - Moreau theorem coherent risk measures can be represented in the dual form as maximum of expected values. Consequently optimization problems involving coherent risk measures can be written in a minimax form. This suggests that such risk averse problems can be formulated as risk neutral problems with respect to an appropriate worst case probability distribution (cf., [15, Remarks 24-25]). Although such worst case probability distribution is not known a priori, we are going to demonstrate that in some cases it can be approximated in a computationally feasible way.

From several points of view a natural example of coherent risk measures is the so-called Average Value at Risk (${\sf AV@R}$) (under different names, such as Conditional Value at Risk, Expected Shortfall, Expected Tail Loss, this was discovered and rediscovered in various equivalent forms by several authors over many years). A nested risk averse approach, using convex combinations of the expectation and ${\sf AV@R}$, was suggested in [16] for controlling high costs in planning of hydropower generation. In that approach risk of high costs is controlled by imposing an appropriate penalization on such high costs at every stage of the decision process conditional on observed realizations of the random data process. One of the criticisms of such risk averse approach is that the corresponding objective is formulated in a nested form and is difficult for an intuitive interpretation.

The contribution of this paper is twofold. We demonstrate that in some situations it is possible to reformulate the considered risk averse problem in a risk neutral form by making an appropriate change of the probability measure. This leads to an intuitive interpretation of controlling the risk by giving higher weights to “bad scenarios”. The idea of constructing scenario trees with extreme (bad) scenarios was considered before (e.g., [17, Chapter 2]). In that respect our approach is quite different. We relate it to the modern risk averse approach to MSP and blend it with the Stochastic Dual Dynamic Programming (SDDP) type algorithm. In particular this allows construction of lower and upper numerical bounds, for the constructed risk neutral problem, following the standard risk neutral methodology of the SDDP method. Also our approach is different from the approach based on extended polyhedral risk measures which allows reformulation of the respective risk-averse problems as risk-neutral problems (with additional constraints and variables) (cf., [7]).

Another contribution of the suggested “change of the probability measure” method is an improvement of the rate of convergence of the straightforward risk averse method. The idea is somewhat related to the classical Importance Sampling techniques although is not exactly the same. The standard Importance Sampling methodology is aimed at reducing variance of the respective sample estimates, while we are concerned with an approximation of the corresponding risk averse problem. An intuitive explanation of our approach is that by generating “bad” scenarios more often one becomes more conservative and risk aware in his/her decisions.

This paper is organized as follows. In the next section we discuss a static case of the risk averse stochastic programming. We show how in some situations the corresponding worst case distribution can be computed. In section 3 we extend this to a multistage setting. In section 4 we discuss a risk averse variant of the Stochastic Dual Dynamic Programming (SDDP) algorithm and its reformulation in a risk neutral form. In section 5 we give a numerical example based on the Brazilian Interconnected Power System operation planning problem. Finally section 6 is devoted to concluding remarks.

2 Static case

Consider the following risk averse stochastic program

[TABLE]

where $(\Omega,{\cal F},P)$ is a probability space, ${\cal X}\subset{\mathbb{R}}^{n}$, $F:{\mathbb{R}}^{n}\times\Omega\to{\mathbb{R}}$, and $\rho:{\cal Z}\to{\mathbb{R}}$ is a coherent risk measure defined on a linear space ${\cal Z}$ of random variables $Z:\Omega\to{\mathbb{R}}$. We assume that for every $x\in{\cal X}$, random variable $F_{x}(\omega)=F(x,\omega)$ belongs to ${\cal Z}$. We also assume that problem (2.1) is convex, i.e., the set ${\cal X}$ is convex and $F(x,\omega)$ is convex in $x$ for a.e. $\omega\in\Omega$. In particular we deal with risk measures of the form

[TABLE]

where ${\cal Z}=L_{1}(\Omega,{\cal F},P)$ paired with its dual space ${\cal Z}^{*}=L_{\infty}(\Omega,{\cal F},P)$, and

[TABLE]

In the above variational form, ${\sf AV@R}$ was defined in [12] under the name “Conditional Value at Risk”.

By the Fenchel - Moreau theorem, real valued coherent risk measure $\rho$ has dual representation (cf., [13])

[TABLE]

where ${\mathfrak{A}}\subset{\cal Z}^{*}$ is a convex weakly∗ compact set of probability density functions and

[TABLE]

is the expectation with respect to the probability measure $dQ=\zeta dP$. Hence problem (2.1) can be written in the following minimax form

[TABLE]

A dual of problem (2.4) is obtained by interchanging the ‘min’ and ‘max’ operators:

[TABLE]

A point $(\bar{x},\bar{\zeta})\in{\cal X}\times{\mathfrak{A}}$ is said to be a saddle point of the above minimax problems if

[TABLE]

Under mild regularity conditions the minimax problem (2.4) has a saddle point $(\bar{x},\bar{\zeta})\in{\cal X}\times{\mathfrak{A}}$. Then $\bar{x}$ is an optimal solution of problem (2.4), $\bar{\zeta}$ is an optimal solution of problem (2.5), optimal values of problems (2.4) and (2.5) are equal to each other and are equal to the optimal value of the following problem

[TABLE]

It follows from the first inequality in (2.6) that if $\bar{x}$ is an optimal solution of problem (2.1) (i.e., $(\bar{x},\bar{\zeta})$ is a saddle point), then $\bar{x}$ is also an optimal solution of problem (2.7). That is, the set of optimal solutions of problem (2.7) contains the set of optimal solutions of problem (2.1) (it can happen that the set of optimal solutions of problem (2.7) is larger than the set of optimal solutions of problem (2.1)).

That is, risk averse problem (2.1) can be formulated as risk neutral problem (2.7) with respect to the “worst” probability measure $dQ=\bar{\zeta}dP$. Of course $\bar{\zeta}$ is not known, its evaluation requires solution of the minimax problem (2.5). Nevertheless this gives us a direction for constructing approximation of problem (2.7). For $Z=F_{\bar{x}}$ we have that

[TABLE]

Recall that for $Z\in{\cal Z}$,

[TABLE]

(cf. [15, eq.(6.43), p.265]).

Consider risk measure (2.2). We have that

[TABLE]

with {\mbox{\boldmath1}}(\cdot)\equiv 1 and (cf. [15, eq.(6.74), p.273])

[TABLE]

where

[TABLE]

Suppose further that the space $\Omega=\{\omega_{1},...,\omega_{N}\}$ is finite equipped with equal probabilities $p_{i}=1/N$, $i=1,...,N$. Then denoting $Z_{i}=Z(\omega_{i})$, random variable $Z:\Omega\to{\mathbb{R}}$ can be identified with vector $(Z_{1},...,Z_{N})\in{\mathbb{R}}^{N}$. In that case

[TABLE]

where $\bar{u}={\sf V@R}_{\alpha}(Z)$. Let $Z_{(1)}\leq...\leq Z_{(N)}$ be values $Z_{i}$, $i=1,...,N$, arranged in the increasing order. Then $\bar{u}=Z_{(\kappa)}$, where $\kappa:=\lceil(1-\alpha)N\rceil$ with $\lceil a\rceil$ denoting the smallest integer $\geq a$. We assume that $N$ is large enough so that $\kappa\leq N-1$. Then we can write

[TABLE]

where

[TABLE]

Note that ${\mathfrak{q}}_{i}\geq 0$, $i=1,...,N$, and $\sum_{i=1}^{N}{\mathfrak{q}}_{i}=1$. That is, in order to find the worst probability measure $Q$ we only need to identify the “bad” outcomes of the distribution of $F_{\bar{x}}$, i.e., which values $Z_{i}=F_{\bar{x}}(\omega_{i})$ are larger than ${\sf V@R}_{\alpha}(F_{\bar{x}})$. We sometimes write ${\mathfrak{q}}_{i}={\mathfrak{q}}_{i}(Z)$ for values defined in (2.14) associated with vector $Z=(Z_{1},...,Z_{N})$.

3 Multistage case

Consider a multistage risk averse stochastic programming problem given in the following nested form (cf., [15])

[TABLE]

driven by the random data process $\xi_{1},\xi_{2},{\dots},\xi_{T}$. Here $x_{t}\in{\mathbb{R}}^{n_{t}}$, $t=1,{\dots},T$, are decision variables, $f_{t}:{\mathbb{R}}^{n_{t}}\times{\mathbb{R}}^{d_{t}}\to{\mathbb{R}}$ are continuous functions, ${\cal X}_{t}:{\mathbb{R}}^{n_{t-1}}\times{\mathbb{R}}^{d_{t}}\rightrightarrows{\mathbb{R}}^{n_{t}}$, $t=2,{\dots},T$, are measurable closed valued multifunctions and $\rho_{t|\xi_{[t-1]}}$ are conditional coherent risk mappings (we use notation $\xi_{[t]}:=(\xi_{1},...,\xi_{t})$ for the history of the data process). As the main example we consider the following conditional counterpart of the risk measure (2.2):

[TABLE]

The first stage data, i.e., the vector $\xi_{1}$, the function $f_{1}:{\mathbb{R}}^{n_{1}}\to{\mathbb{R}}$, and the set ${\cal X}_{1}\subset{\mathbb{R}}^{n_{1}}$ are deterministic.

We assume that problem (3.1) is convex, i.e., functions $f_{t}(\cdot,\xi_{t})$ and sets ${\cal X}_{t}(x_{t-1},\xi_{t})$ are convex. It is said that the multistage problem (3.1) is linear if the objective functions and the constraint functions are linear, that is

[TABLE]

Here, $\xi_{1}:=(c_{1},A_{1},b_{1})$ is known at the first stage (and hence is nonrandom), and $\xi_{t}:=(c_{t},B_{t},A_{t},b_{t})\in{\mathbb{R}}^{d_{t}}$, $t=2,{\dots},T$, are data vectors some (or all) elements of which can be random. Linear problems are convex.

Problem (3.1) can be written in the following equivalent form

[TABLE]

where $\Pi$ is the set of policies $\pi=\left(x_{1},x_{2}(\xi_{[2]}),...,x_{T}(\xi_{[T]})\right)$ satisfying the feasibility constraints of problem (3.1), and $\bar{\rho}$ is the composite risk measure (cf., [15, p.318])

[TABLE]

The composite risk measure $\bar{\rho}$ is given in the nested form (3.5) and could be quite complicated. As it was already mentioned in the introduction, this may raise an objection of an intuitive interpretation of the overall objective of the risk averse formulation (3.4).

Anyway risk measure $\bar{\rho}$ is coherent and has the dual representation

[TABLE]

where ${\mathfrak{M}}$ is a set of probability measures (distributions) of $\xi_{[T]}$ absolutely continuous with respect to the reference probability measure of the data process. Therefore problem (3.4) can be written in the following minimax form

[TABLE]

The dual of problem (3.7) is the problem

[TABLE]

Under mild regularity conditions optimal values of problems (3.7) and (3.8) are equal to each other. Suppose further that problem (3.8) has an optimal solution $\bar{Q}$. Then problem (3.1) has the same optimal value as the risk neutral problem

[TABLE]

and the set of optimal solutions of problem (3.1) is contained in the set of optimal solutions of problem (3.9). When the number of scenarios is finite, i.e., the data process can be represented by a finite scenario tree, the change of measure can be described in a constructive way (cf., [15, Remarks 24-25, pp.314-315]).

We assume in the remainder of this section that the data process $\xi_{t}$ is stagewise independent, i.e., random vector $\xi_{t+1}$ is independent of $\xi_{[t]}$, $t=1,...,T-1$ (although $\xi_{1}$ is supposed to be deterministic, we include it for uniformity of the notation.). Suppose further that the conditional risk mappings $\rho_{t|\xi_{[t-1]}}$ are given as the conditional counterparts of coherent risk measures $\rho_{t}$. In the stagewise independent case the joint probability distribution of $(\xi_{1},...,\xi_{T})$ is determined by the marginal distributions of each $\xi_{t}$, $t=1,...,T$, and the corresponding cost-to-go (value) functions can be written as

[TABLE]

where

[TABLE]

$t=1,...,T$, with ${\cal V}_{T+1}(\cdot)\equiv 0$.

Let

[TABLE]

Note that $\bar{x}_{t}$, $t=2,...,T$, is a function of $\bar{x}_{t-1}$ and $\xi_{t}$, and that the policy $(\bar{x}_{1},...,\bar{x}_{T})$ is an optimal solution of the corresponding multistage problem (3.1). For each $\rho_{t+1}$ we can consider its dual representation of the form (2.3)

[TABLE]

with the corresponding set ${\mathfrak{A}}_{t+1}$ of density functions. Consider a saddle point of the minimax problem

[TABLE]

For $x_{t}=\bar{x}_{t}$ we need to solve the problem

[TABLE]

in order to find the component $\bar{\zeta}_{t+1}$ of the corresponding saddle point $(\bar{x}_{t},\bar{\zeta}_{t+1})$.

Suppose now that marginal distribution of $\xi_{t}$, $t=2,...,T$, is discretized by generating $N$ points $\xi_{t}^{1},...,\xi_{t}^{N}$, each assigned with the same probability $1/N$. Suppose further that $\rho_{t}=\rho$, for all $t$, with $\rho$ given in the form (2.2). Then as it was discussed in section 2, in order to find the corresponding worst case density (worst case distribution) we only need to identify the “bad” outcomes of the distribution of $V_{t+1}(\bar{x}_{t},\xi_{t+1})$ (see (2.14)), i.e., which values $V_{t+1}(\bar{x}_{t},\xi^{j}_{t+1})$ are larger than ${\sf V@R}_{\alpha}(V_{t+1}(\bar{x}_{t},\xi_{t+1}))$. Recall that $\bar{x}_{t}$, $t=2,...,T$, is a function of $\bar{x}_{t-1}$ and $\xi_{t}$. Nevertheless if ${\sf V@R}_{\alpha}(V_{t+1}(\bar{x}_{t},\xi_{t+1}))$ is stable with respect to different realizations of $\bar{x}_{t}$, then it would be possible to assign weights (probabilities) to $\xi_{t}^{j}$, at every stage $t$ of the decision process, in such a way that the constructed risk neutral problem will be equivalent to the original risk averse multistage problem. We will discuss such an example in the following sections.

4 Stochastic Dual Dynamic Programming algorithm

In the risk neutral setting the Stochastic Dual Dynamic Programming (SDDP) algorithm was introduced in Pereira and Pinto [8]. Its origins can be traced to the nested decomposition algorithm of Birge [2]. A risk averse variant of the SDDP method, based on a convex combination of the expectation and ${\sf AV@R}$, was introduced in [14] and implemented in [10]. For extended polyhedral risk measures a variant of the SDDP method is discussed in [7]. Convergence of the SDDP type algorithms is discussed in [4],[5],[6], and references therein. For a view of the SDDP method as a form of approximate dynamic programming we refer to [powell].

Consider the multistage stochastic programming problem (3.1). We assume that the problem is linear, i.e., the data are given in the form (3.3). Moreover we assume that the data process $\xi_{t}$ is stagewise independent and has a finite number of scenarios. That is, the marginal distribution of $\xi_{t}$, $t=2,...,T$, has $N$ realizations $\xi_{t}^{1},...,\xi_{t}^{N}$ each having the same probability $1/N$ (for the sake of simplicity we assume that the number $N$ of discretization points at each stage $t$ is the same). The SDDP algorithm is a cutting plane type method designed to solve such multistage convex stochastic programming problems. Being an iterative approach, the SDDP algorithm progressively refines lower piecewise linear approximations of ${\cal V}_{t+1}(\cdot)$. It has two major steps at each iteration.

Forward Step

At iteration $m$, suppose for each $t=1,\ldots,T-1$, we have a finite set $\mathcal{S}_{t+1}^{m}$ of affine minorants $s(\cdot)$ of ${\cal V}_{t+1}(\cdot)$, then $\mathfrak{V}_{t+1}^{m}(\cdot):=\max_{s\in\mathcal{S}_{t+1}^{m}}s(\cdot)$ is a lower piecewise linear approximation of ${\cal V}_{t+1}(\cdot)$, and

[TABLE]

is a lower convex approximation of $V_{t}(\cdot,\xi_{t})$ for each $\xi_{t}$.

In the forward step of the algorithm a random realization $\xi_{t}^{j_{t}}$, $j_{t}\in\{1,....,N\}$, from the distribution of $\xi_{t}$, $t=2,...,T$, is sampled (we refer to this as the uniform sampling). In the case that each realization $\xi_{t}^{j}$ of $\xi_{t}$ has the same probability $1/N$, each scenario (sample path) $\{\xi_{t}^{j_{t}}\}_{t=2,\ldots,T}$ of the data process is sampled with probability $1/N^{T-1}$. At stage $t$ the corresponding trial point $x_{t}^{m}$ is defined to be the optimal solution

[TABLE]

Note that $f_{1}(x_{1}^{m})+\mathfrak{V}_{2}^{m}(x_{1}^{m})$ gives a lower bound of the total cost.

Remark 4.1

When the problem is risk neutral (i.e., $\lambda=0$), the measure $\bar{Q}$ in (3.9) is given by a discrete measure that assigns probability $1/N^{T-1}$ to each scenario $\{\xi_{t}^{j_{t}}\}$. Viewing $\{x_{t}^{m}\}$ as functions of $\{\xi_{t}^{j_{t}}\}$, the quantity $\mathbb{E}_{\bar{Q}}[\sum_{t=1}^{T}f_{t}(x_{t}^{m}(\xi_{t}),\xi_{t})]$ is then an upper bound of the (optimal) total cost, and $\sum_{t=1}^{T}f_{t}(x_{t}^{m},\xi_{t}^{j_{t}})$ is the corresponding unbiased estimator. By generating a number of scenarios one can construct a confidence interval for the upper bound. In the risk averse setting construction of the corresponding upper bound is more involved. A straightforward analogue of the statistical upper bound, used in the risk neutral case, does not work in the risk averse setting (cf., [3]). An approach to constructing (nonstochastic) upper bounds, which is also applicable in the risk averse setting, was suggested in [11]. **

Backward Step

Given trial points $\{x_{t}^{m}\}_{t=1,\ldots,T-1}$ at iteration $m$, the approximations $\mathfrak{V}_{t}^{m}(\cdot)$ are refined sequentially from $t=T$ down to $t=2$. By equation (2.13), we have

[TABLE]

where $\underline{V}_{t}^{m}(x_{t-1},\xi_{t}^{(1)})\leq\ldots\leq\underline{V}_{t}^{m}(x_{t-1},\xi_{t}^{(N)})$ are values $\underline{V}_{t}^{m}(x_{t-1},\xi_{t}^{j}))$, $j=1,\ldots,N$, arranged in the increasing order, and $\kappa:=\lceil(1-\alpha)N\rceil$. The corresponding affine minorant $s_{t}^{m}(\cdot)$ of $\rho_{t}(\underline{V}_{t}^{m+1}(\cdot,\xi_{t}))$ at $x_{t-1}^{m}$ is constructed by computing a subgradient $g_{t}^{j}$ of $\underline{V}_{t}^{m+1}(\cdot,\xi_{t}^{j})$ at $x_{t-1}^{m}$ (cf., [16, Section 4.1]). Observe that $s_{t}^{m}(\cdot)$ is also an affine minorant of $\mathcal{V}_{t}(\cdot)$, since $\rho_{t}(\underline{V}_{t}^{m+1}(\cdot,\xi_{t}))$ is a lower convex approximation of $\mathcal{V}_{t}(\cdot)$. Next, define $\mathcal{S}_{t}^{m+1}:=\mathcal{S}_{t}^{m}\cup\{s_{t}^{m}\}$. Then $\mathfrak{V}_{t}^{m+1}(\cdot)=\max\{\mathfrak{V}_{t}^{m}(\cdot),s_{t}^{m}(\cdot)\}$ is a refined lower piecewise linear approximation of $\mathcal{V}_{t}(\cdot)$.

We now discuss how to generate the subgradients $g_{t}^{j}$ when the multistage problem is linear (see (3.3)). Suppose $\xi_{t}^{j}=(c_{t}^{j},B_{t}^{j},A_{t}^{j},b_{t}^{j})$, and $\mathfrak{V}_{t+1}^{m+1}(\cdot)$ has been generated, where $\mathfrak{V}_{T+1}^{m+1}(\cdot)\equiv 0$. By definition (4.16),

[TABLE]

A subgradient of $\underline{V}_{t}^{m+1}(\cdot,\xi_{t}^{j})$ at $x_{t-1}^{m}$ is given by $g_{t}^{j}=-(B_{t}^{j})^{\top}\pi_{t}^{j}$, where $\pi_{t}^{j}$ is an optimal solution to the dual problem of

[TABLE]

Problem (4.20) can be reformulated as the linear programming problem

[TABLE]

The dual problem of (4.21) is also a linear programming problem, hence a dual optimal solution $\tilde{\pi}_{t}^{j}$ can be obtained efficiently, and $\pi_{t}^{j}$ is exactly the entries of $\tilde{\pi}_{t}^{j}$ corresponding to the constraints $B_{t}^{j}x_{t-1}^{m}+A_{t}^{j}x_{t}=b_{t}^{j}$.

Remark 4.2

The backward steps can be carried out with arbitrary (feasible) points $x_{t}^{m}$, $t=1,\ldots,T-1$. The trial points generated at the forward steps are one of the possibilities. As we shall see later, there could be more efficient ways to generate the trial points. **

Remark 4.3

Depending on whether $\lambda>0$ or $\lambda=0$ we refer to the SDDP algorithm as risk averse or risk neutral, respectively.

4.1 Identifying “bad” outcomes

As mentioned in the last paragraph of section 3, we can find the corresponding worst case probability density by identifying the bad outcomes. To identify the bad outcomes, we need to verify the order of values of $V_{t}(\bar{x}_{t-1},\xi_{t}^{j})$, $j=1,\ldots,N$, at each stage $t=2,\ldots,T$. Although the required values $V_{t}(\bar{x}_{t-1},\xi_{t}^{j})$ are not available, we have access to the approximations $\underline{V}_{t}^{m+1}(x_{t-1}^{m},\xi_{t}^{j})$ at the $m$-th iteration of the SDDP algorithm. If the approximations are good enough, by continuity argument the order of $\underline{V}_{t}^{m+1}(x_{t-1}^{m},\xi_{t}^{j})$ is more or less the same as the order of $V_{t}(\bar{x}_{t-1},\xi_{t}^{j})$, $j=1,...,N$. Since the approximations are improved in each iteration of the SDDP algorithm, we expect the order of $\underline{V}_{t}^{m+1}(x_{t-1}^{m},\xi_{t}^{j})$ to stabilize after a certain number of iterations. In particular, if the SDDP algorithm is run for a sufficient number of iterations, then those bad outcomes should be the outcomes that correspond to higher values of $\underline{V}_{t}^{m+1}(x_{t-1}^{m},\xi_{t}^{j})$ (in each iteration) the most frequently.

Following this idea, we try to identify the bad outcomes via frequency, i.e., by counting the number of iterations that values

[TABLE]

are large at a given stage $t$. Formally, let $\nu_{t,(1)}^{m+1}\leq\ldots\leq\nu_{t,(N)}^{m+1}$ be these values in the increasing order, $\kappa:=\lceil(1-\alpha)N\rceil$ and

[TABLE]

be the corresponding quantile. Let $|A|$ denote the cardinality of a set $A$, and define

[TABLE]

In plain words, for each $t$ and $j_{0}$, $W_{t}^{j_{0}}$ is the frequency (i.e., the number of iterations) that $\underline{V}_{t}^{m+1}(x_{t-1}^{m},\xi_{t}^{j_{0}})$ is one of the $\lceil\alpha N\rceil$-largest values in $\{\underline{V}_{t}^{m+1}(x_{t-1}^{m},\xi_{t}^{j})\}_{1\leq j\leq N}$. Ideally, the bad outcomes $\xi_{t}^{j}$ should correspond to a high value $W_{t}^{j}$. For the operation planning problem studied in section 5, numerical experiments indicate that $W_{t}^{j}$ can be as high as $2900$, after a total of $3000$ iterations, in most of the stages.

We next discuss how to assign probability weights to outcomes $\xi_{t}^{j}$ based on values of $W_{t}^{j}$. Let $W_{t}^{(1)}\leq\ldots\leq W_{t}^{(N)}$ be values $W_{t}^{j}$ arranged in the increasing order and $q_{t}^{j}:={\mathfrak{q}}_{j}(W_{t})$, $j=1,\ldots,N$, be values defined in (2.14) associated with vector $W_{t}=(W^{1}_{t},\ldots,W_{t}^{N})$, then we assign weights $q_{t}^{j}$ to $\xi_{t}^{j}$. If $W_{t}^{(j)}$ are not in the strictly increasing order (i.e., some of these values are equal to each other), we still may sort $W_{t}^{j}$ to ensure that $\sum_{j=1}^{N}q_{t}^{j}=1$ for $t=1,\ldots,T$.

Remark 4.4

Suppose the frequencies $W_{t}^{j}$ are given, and the weights $q_{t}^{j}$ are assigned to outcomes $\xi_{t}^{j}$ based on the values of $W_{t}^{j}$. Let each scenario $\{\xi_{t}^{j_{t}}\}_{t=2,\ldots,T-1}$ be sampled at forward steps with probability $\prod_{t=2}^{T-1}q_{t}^{j_{t}}$ (we refer to this as the biased sampling). Then the forward steps produce trial points $\{x_{t}^{m}\}_{t=1,\ldots,T-1}$ that are different from the ones obtained under uniform sampling (as described in (4.17)). It turns out that the convergence of the risk averse SDDP algorithm is improved when we utilize those points at backward steps.

The improvement in convergence can be explained in the following way. By (2.13) we have that

[TABLE]

with $q_{t,j}^{m}={\mathfrak{q}}_{j}(\nu_{t}^{m})$. When the order of $\{\underline{V}_{t}^{m}(x_{t-1},\xi_{t}^{j})\}$ is stable, sampling $\xi_{t}^{j}$ that corresponds to higher weights $q_{t,j}$ increases chances of appearing bad scenarios. This is reminiscent of the Importance Sampling (IS) techniques, although is not exactly the same. **

Remark 4.5** (Change-of-measure risk neutral problem)**

We can construct a new risk neutral problem, referred to as the change-of-measure risk neutral problem, such that the outcomes $\{\xi_{t}^{j}\}$ at stage $t$ are assigned the respective weights $\{q_{t}^{j}\}$. Such problem can be solved by applying the risk neutral SDDP algorithm with uniform sampling at forward steps. Here, uniform sampling means that each scenario $\{\xi_{t}^{j_{t}}\}_{t=2,\ldots,T-1}$ is sampled with probability $\prod_{t=2}^{T-1}q_{t}^{j_{t}}$ rather than $1/N^{T-2}$, since the weights of $\xi_{t}^{j}$ are $q_{t}^{j}$ instead of $1/N$.

In section 5, we apply the risk neutral SDDP algorithm (with $\lambda=0$) to the change-of-measure risk neutral problem constructed for the operation planning problem. The risk neutral SDDP algorithm yields similar approximations as the ones obtained by applying the risk averse SDDP algorithm to the original risk averse formulation of the operation planning problem. This offers an explanation of the nested risk averse formulation of the objective function, that such formulation is equivalent to the risk neutral problem where the risk is implicitly controlled by assigning higher weights to “bad scenarios”. **

4.2 Dynamic Biased Sampling

In Remark 4.4, we have discussed how a biased sampling method for generating scenarios (thus trial points) at forward steps can help to improve convergence of the risk averse version of the SDDP algorithm. However, those bad outcomes are not known a priori. In order to acquire the frequencies $W_{t}^{j}$, we need to run the risk averse SDDP algorithm with uniform sampling first, which can be time consuming given the size of the problem. For the purpose of improving the convergence of the risk averse SDDP algorithm, we can instead employ a dynamic biased sampling method (for generating scenarios at forward steps) that allocates more weights to a dynamic set of outcomes. Such set gets updated in each iteration throughout the algorithm. As illustrated in section 5, the performance of the risk averse SDDP algorithm equipped with dynamic biased sampling is similar to the performance of the algorithm with biased sampling.

We next discuss how the dynamic set is updated. Consider variables $\gamma_{t}^{j},t=1,\ldots,T,j=1,\ldots,N$, that represents the adjusted frequency of outcomes $\xi_{t}^{j}$. Initially, set $\gamma_{t}^{j}=0$ for all $t,j$. At the backward step of iteration $m$, we increase $\gamma_{t}^{j}$ by $1$ if $\underline{V}_{t}^{m+1}(x_{t-1}^{m},\xi_{t}^{j})\geq\nu_{t,(\kappa)}^{m+1}$. We then set $\gamma_{t}^{j}\leftarrow\gamma_{t}^{j}\cdot\frac{m}{m+1}$ for all $t$ and $j$, and proceed to the next iteration. The adjustment weakens the impact of early iterations, where the approximations are still inaccurate. The dynamic set of bad outcomes in each stage $t$ are the outcomes with the highest adjusted frequency so far, and weights are assigned according to (2.14). When the bad outcomes are stable, the dynamic set should coincide with the set of those outcomes. A complete description of the algorithm can be found in Appendix A.

5 Numerical Experiments

The numerical experiments are performed on an aggregated representation of the Brazilian Interconnected Power System operation planning problem with historical data as of January 2011. The study horizon is 60 stages and the total number of considered stages is $T=120$. The scenario tree has 100 realizations at every stage (i.e., $N=100$), and each realization $\xi_{t}^{j},t=2,\ldots,120,j=1,\ldots,100$, has the same weight $1/N=1/100$. The random data process is represented by four dimensional vectors of monthly inflows, aggregated in four regions, and modeled as the first order periodical time series process. The total number of state variables is 8. We refer to [16] for the detail description of the model (see also Appendix B). We implement the risk averse SDDP algorithms with different sampling methods (for generating scenarios and thus trial points). We also implement the risk neutral SDDP algorithm for the constructed change-of-measure risk neutral problem (see Remark 4.5). Both implementations were written in C++ and using Gurobi 8.1. Dual simplex was used as a default method for the LP solver.

We conduct the numerical experiments with parameters $\lambda=0.2$ and $\alpha=0.05$ (see (3.2)) in the following three steps:

The first step is to run the risk averse SDDP algorithm with uniform sampling, referred to as “raus” (where each scenario $\{\xi_{t}^{j_{t}}\}_{t=2,\ldots,T-1}$ is sampled with probability $1/N^{T-2}$ at forward steps) and to identify the bad outcomes as discussed in subsection 4.1. 2. 2.

The second step is to run the risk averse SDDP algorithm with biased sampling, referred to as “rabs” (where each scenario $\{\xi_{t}^{j_{t}}\}_{t=2,\ldots,T-1}$ is sampled with probability $\prod_{t=2}^{T-1}q_{t}^{j_{t}}$, see Remark 4.4) and with dynamic biased sampling “radbs” (see section 4.2). We then compare lower bounds of the optimal objective value of the risk averse problem produced by the SDDP algorithm with different sampling methods (recall that the lower bound is given by the quantity $f_{1}(x_{1}^{m})+\mathfrak{V}_{2}^{m}(x_{1}^{m})$). 3. 3.

The final step is to solve the change-of-measure risk neutral problem by the risk neutral SDDP algorithm with uniform sampling, referred to as “nrn” (where each scenario $\{\xi_{t}^{j_{t}}\}_{t=2,\ldots,T-1}$ is sampled with probability $\prod_{t=2}^{T-1}q_{t}^{j_{t}}$, see Remark 4.5). In addition, we compare the policies corresponding to the risk averse problem and the change-of-measure risk neutral problem.

In addition, we perform numerical experiments in the extreme case where $\lambda=0.5$ and $\alpha=0.05$ to examine the effectiveness of dynamic biased sampling in improving the convergence of the lower bounds.

Experiment with $\lambda=0.2,\alpha=0.05$.

We first run the risk averse SDDP algorithm with uniform sampling for 3000 iterations. As discussed in subsection 4.1, the bad outcomes are identified via frequencies $W_{t}^{j}$, and the weights $q_{t}^{j}$ are assigned accordingly. Here, $\kappa=\lceil(1-\alpha)N\rceil=95$, and more weights are assigned to the outcomes corresponding to the top five frequencies.

Figure 1 plots the top five frequencies at each stage after 3000 iterations. In particular, the $i$-th curve corresponds to the $i$-th highest frequency by stage, which are the points $\{(t,W_{t}^{(N-i+1)})\}_{t=2,\ldots,T}$ in the notation in subsection 4.1. A high frequency suggests the corresponding outcome is indeed a bad outcome and should be assigned more weight in the worst case probability density. Observe that the frequency of the top three curves are almost the same as the total number of iterations. The fluctuations of frequency of the $4$ and $5$-th curve were partly due to inaccuracy of the approximations in the early iterations, and the corresponding outcomes become stable as the number of iterations increases. Such evidence supports the validity of the change-of-measure approach.

Figure 2 contains two plots of lower bounds of the optimal objective value of the considered risk averse problem produced by the risk averse SDDP algorithm with different sampling methods (for generating scenarios thus trial points at forward steps), i.e., for each sampling method, we plot the points $(m,f_{1}(x_{1}^{m})+\mathfrak{V}_{2}^{m}(x_{1}^{m}))$ by iteration $m$ (see also section 4).

The plot on the left of Figure 2 compares “raus”, “rabs”, and “radbs”, which correspond to the risk averse SDDP algorithm where the trial points are generated under uniform, biased, and dynamic biased sampling, respectively. It shows the convergence of “rabs” and “radbs” are almost the same, which are much better than that of “raus”. As discussed in subsection 4.2, “rabs” requires weights obtained by first running “raus”, which can be time consuming given the size of the problem. This makes “radbs” a better choice for the purpose of improving the convergence of the risk averse SDDP algorithm. Nevertheless, the approximations generated by “rabs” and “radbs” are similar, and they are better than the one generated by “raus”. In the remaining section, we use the approximation $\{\mathfrak{V}_{t}^{A}\}$ generated by “radbs” after 3000 iterations as the reference approximation for the risk averse problem.

The plot on the right of Figure 2 compares “raus”, “radbs ”, and “raus+bs” in varying iterations, where “raus+bs” is the risk averse SDDP algorithm equipped with uniform sampling in the first 3000 iterations and switched to biased sampling afterwards. We see that the lower bound generated by “raus” in 8000 iterations is still much worse than the one generated by “rabs” in 3000 iterations, and the gap closes slowly. Besides, “raus+bs” takes about 6000 iterations to reach what “radbs” achieves in 3000 iterations.

As the final step, we solve the change-of-measure risk neutral problem by “nrn” (see step 3). Let $\{\mathfrak{V}_{t}^{N}\}$ denote the approximations produced by “nrn” after 3000 iterations. We next compare the policies generated by the approximations $\{\mathfrak{V}_{t}^{A}\}$ and $\{\mathfrak{V}_{t}^{N}\}$. Formally, for a given scenario $\{\xi_{t}^{j_{t}}\}_{t=1,\ldots,T}$ and approximations $\{\mathfrak{V}_{t}\}_{t=1,\ldots,T}$, a policy generated by the approximations is $\{x_{t}(\xi_{t}^{j_{t}})\}_{t=1,\ldots,T}$ such that

[TABLE]

The policies are the decisions to be implemented given the scenario. If $\{\mathfrak{V}_{t}^{A}\}$ and $\{\mathfrak{V}_{t}^{N}\}$ yield similar policies for each scenario, then the risk averse formulation is similar to the change-of-measure risk neutral formulation, at least from a numerical perspective. Instead of comparing policies for all scenarios, we randomly sample 3000 scenarios and compare the policies by plotting their paths (or distribution).

For operation planning problem, the set of stage variables is the same across all stages, and each stage variable at stage $t$ corresponds to an entry of $x_{t}(\xi_{t}^{j_{t}})$. We sample 3000 scenarios uniformly at random and generate a fanplot for paths of each variable and of stage objective value (explained below). In Figures 3, 4 and 5, we present fanplots that exhibit typical behavior of the variables. The dark area on the fanplots is where the paths of the variables are highly concentrated, whereas the light area is the opposite. The orange curve on the plots is the average of the variables by stage.

The plots on the left of Figures 3, 4 and 5 correspond to the risk averse problem, whereas the ones on the right correspond to the constructed risk neutral problem. We can see from the plots that concentrations of policies of the same variable generated by $\{\mathfrak{V}_{t}^{A}\}$ and $\{\mathfrak{V}_{t}^{N}\}$ are similar.

Figure 6 plots concentrations of the stage objective value $f_{t}(x_{t}(\xi_{t}^{j_{t}}),\xi_{t}^{j_{t}})$ for policies $x_{t}(\xi_{t}^{j_{t}})$ generated by $\{\mathfrak{V}_{t}^{A}\}$ and $\{\mathfrak{V}_{t}^{N}\}$, respectively. Figures 3 to 6 indicate that the risk averse and the constructed risk neutral problem yield similar policies and objective values, hence the risk averse formulation and the change-of-measure risk neutral formulation are similar.

Figure 7 plots the lower bounds and statistical upper bounds of the optimal objective value of the change-of-measure risk neutral problem. Here, the statistical upper bound at iteration $m$ is defined to be the quantity $\frac{1}{m}\sum_{i=1}^{m}\sum_{t=1}^{T}f_{t}(x_{t}^{i}(\xi_{t}^{j_{t}^{i}}),\xi_{t}^{j_{t}^{i}})$, where $\xi_{t}^{j_{t}^{i}}$, $t=1,\ldots,T$, are the scenarios generated at iteration $i$ (note that it is different from the statistical upper bound introduced in [8]). After 3000 iterations the gap between the lower and upper bounds is about 14%. One advantage of applying the SDDP algorithm to the constructed risk neutral problem is the possibility of incorporating the termination criterion based on the lower and upper bounds.

Experiment with $\lambda=0.5,\alpha=0.05$.

The case $\lambda=0.5$ leads to an extreme risk aversion. The parameter $\lambda=0.5$ was chosen instead of some larger values (e.g. $\lambda=1$) to better reflect the reality.

Figure 8 compares the convergence of lower bounds generated by “raus”, “radbs”, and two variants of “radbs”, namely “radbsm1” and “radbsm2”. Those two variants differ from “radbs” in the adjustment $\gamma_{t}^{j}\leftarrow\gamma_{t}^{j}\cdot\frac{m}{m+1}$ in iteration $m$; “radbsm1” simply does not have such adjustment, whereas the adjustment in “radbsm2” is replaced by $\gamma_{t}^{j}\leftarrow\gamma_{t}^{j}\cdot(1-0.5^{m})$. From Figure 8, we can see the convergence of the three dynamic biased sampling methods is more or less the same. However, they achieve a striking speed up over “raus”, where their lower bounds are almost twice the lower bound of “raus” after 3000 iterations. Hence, the dynamic biased sampling method gets more effective when $\lambda$ increases. Note in the other extreme case where $\lambda=0$, biased sampling reduces to uniform sampling, thus the dynamic biased sampling approach does not offer any speed up.

6 Conclusions

We discuss the risk averse (multistage) stochastic programming and the idea of “bad” outcomes. We show how to identify the “bad” outcomes and use this with the SDDP method. Numerical experiments were conducted on the Brazilian Interconnected Power System operation planning problem, with the results presented and discussed in section 5. The convergence of lower bounds generated by the risk averse version of the SDDP algorithm with different sampling methods was examined. It was observed that the (dynamic) biased sampling approach for generating trial points considerably improved performance of the lower bounds. We also compared the solutions generated by the risk averse version of the SDDP method and the changed-of-measure risk neutral SDDP method and concluded that the policies and the objective values generated by these two approaches are similar.

Acknowledgement

The authors are indebted to Filipe Goulart Cabral and anonymous referees for constructive comments which helped to improve presentation of the manuscript.

7 Appendix A

8 Appendix B

Long Term Operation Planning Problem

The dynamic programming equation for the long term operation planning problem can be written as

[TABLE]

for $t=1\ldots,T$, where

[TABLE]

and

[TABLE]

with $\lambda\in[0,1]$ and $\alpha\in[0,1]$ being chosen parameters.

The objective function

[TABLE]

represents the sum of the total cost for thermal generation and deficit with $\mathcal{Q}_{t+1}([v_{t+1},a_{[t-p+1,t]}])$, where

• •

$\beta$ is a discount factor;

• •

$K$ is a subsystem set;

• •

$T_{k}$ is the thermal set for subsystem $k$;

• •

$U_{k}$ is the deficit set for subsystem $k$.

The energy balance equation for each reservoir $k$ is

[TABLE]

where

• •

$v_{t}$ is the stored energy in the reservoir at the beginning of stage $t$;

• •

$a_{t}$ is the energy inflow during stage $t$;

• •

$q_{t}$ is the generated energy during stage $t$;

• •

$s_{t}$ is the spilled energy during stage $t$.

The time-series model for the energy inflow is

[TABLE]

where

• •

$a_{t-\nu}$ is the energy inflow during stage $t-\nu$;

• •

$\Phi_{t,\nu}$ is the coefficient of PVAR vector time-series model $t$;

• •

$\eta_{t}$ is the multiplicative noise of PVAR, which independent for each stage $t$.

The load balance equation, in MW month, for each subsystem $k$ and stage $t$ is

[TABLE]

where

• •

$d_{tk}$ is load;

• •

$q_{tk}$ is hydro generation;

• •

$\sum_{j\in T_{k}}g_{tj}$ is thermal generation;

• •

$\sum_{i\in U_{k}}Def_{ti}$ is deficit generation;

• •

$\sum_{\l\in\Omega_{k}}(f_{tlk}-f_{tkl})$ is net energy interchange;

• •

$f_{tlk}$ is the energy flow from subsystem $l$ to subsystem $k$;

• •

$\Omega_{k}$ is the subsystems directly connected to subsystem $k$.

The bounds on variables are

• •

$0\leq v_{t+1}\leq\overline{v}$ is the lower and upper bounds on stored energy;

• •

$0\leq q_{t}\leq\overline{q}$ is the lower and upper bounds on generated energy;

• •

$0\leq s_{t}$ is non-negativity constraint of spilled energy;

• •

$\underline{g}\leq g_{t}\leq\overline{g}$ is the lower and upper bounds on thermal generation;

• •

$0\leq Def_{t}\leq\overline{Def}$ is the lower and upper bounds on energy deficit;

• •

$\underline{f}\leq f_{t}\leq\overline{f}$ is the lower and upper bounds on energy flow.

The main idea of the deficit is to penalize the load cut by a convex piecewise linear cost function which is dependent on the load cut depth. Regarding this approach, it is important to emphasize that the deficit variable with highest associated cost is unbounded above.

For each stage $t$ the decision vector is $x_{t}=(v_{t+1},q_{t},s_{t},g_{t},Def_{t},f_{t},a_{t})$. In the long term operation planning problem the only considered uncertainty is the independent multiplicative noise, that is, $\xi=\eta_{t}$, and the cost-to-go function depends only on $[v_{t},a_{[t-p,t-1]}]$ of last $p$ previous decision $x_{[t-p,t-1]}$. In this way, it is usual to write $Q_{t}([v_{t},a_{[t-p,t-1]}],\eta_{t})$ instead of $Q_{t}(x_{[t-p,t-1]},\eta_{t})$.