Stochastic programming

The formulation of the above two-stage problem assumes that the second-stage data ${\displaystyle \xi }$ is modeled as a random vector with a known probability distribution. This would be justified in many situations. For example, the distribution of ${\displaystyle \xi }$ could be inferred from historical data if one assumes that the distribution does not significantly change over the considered period of time. Also, the empirical distribution of the sample could be used as an approximation to the distribution of the future values of ${\displaystyle \xi }$. If one has a prior model for ${\displaystyle \xi }$, one could obtain a posteriori distribution by a Bayesian update.

To solve the two-stage stochastic problem numerically, one often needs to assume that the random vector ${\displaystyle \xi }$ has a finite number of possible realizations, called scenarios, say ${\displaystyle \xi _{1},\dots ,\xi _{K}}$, with respective probability masses ${\displaystyle p_{1},\dots ,p_{K}}$. Then the expectation in the first-stage problem's objective function can be written as the summation:

$${\displaystyle E[Q(x,\xi )]=\sum \limits _{k=1}^{K}p_{k}Q(x,\xi _{k})}$$

and, moreover, the two-stage problem can be formulated as one large linear programming problem (this is called the deterministic equivalent of the original problem, see section § Deterministic equivalent of a stochastic problem).

When ${\displaystyle \xi }$ has an infinite (or very large) number of possible realizations the standard approach is then to represent this distribution by scenarios. This approach raises three questions, namely:

1. How to construct scenarios, see § Scenario construction;
2. How to solve the deterministic equivalent. Optimizers such as CPLEX, and GLPK can solve large linear/nonlinear problems. The NEOS Server,[6] hosted at the University of Wisconsin, Madison, allows free access to many modern solvers. The structure of a deterministic equivalent is particularly amenable to apply decomposition methods,[7] such as Benders' decomposition or scenario decomposition;
3. How to measure quality of the obtained solution with respect to the "true" optimum.

These questions are not independent. For example, the number of scenarios constructed will affect both the tractability of the deterministic equivalent and the quality of the obtained solutions.

With a finite number of scenarios, two-stage stochastic linear programs can be modelled as large linear programming problems. This formulation is often called the deterministic equivalent linear program, or abbreviated to deterministic equivalent. (Strictly speaking a deterministic equivalent is any mathematical program that can be used to compute the optimal first-stage decision, so these will exist for continuous probability distributions as well, when one can represent the second-stage cost in some closed form.) For example, to form the deterministic equivalent to the above stochastic linear program, we assign a probability ${\displaystyle p_{k}}$ to each scenario ${\displaystyle k=1,\dots ,K}$. Then we can minimize the expected value of the objective, subject to the constraints from all scenarios:

${\displaystyle {\begin{array}{lccccccccccccc}{\text{Minimize}}&f^{\top }x&+&g^{\top }y&+&p_{1}h_{1}^{\top }z_{1}&+&p_{2}h_{2}^{T}z_{2}&+&\cdots &+&p_{K}h_{K}^{\top }z_{K}&&\\{\text{subject to}}&Tx&+&Uy&&&&&&&&&=&r\\&&&V_{1}y&+&W_{1}z_{1}&&&&&&&=&s_{1}\\&&&V_{2}y&&&+&W_{2}z_{2}&&&&&=&s_{2}\\&&&\vdots &&&&&&\ddots &&&&\vdots \\&&&V_{K}y&&&&&&&+&W_{K}z_{K}&=&s_{K}\\&x&,&y&,&z_{1}&,&z_{2}&,&\ldots &,&z_{K}&\geq &0\\\end{array}}}$

We have a different vector ${\displaystyle z_{k}}$ of later-period variables for each scenario ${\displaystyle k}$. The first-period variables ${\displaystyle x}$ and ${\displaystyle y}$ are the same in every scenario, however, because we must make a decision for the first period before we know which scenario will be realized. As a result, the constraints involving just ${\displaystyle x}$ and ${\displaystyle y}$ need only be specified once, while the remaining constraints must be given separately for each scenario.

A common approach to reduce the scenario set to a manageable size is by using Monte Carlo simulation. Suppose the total number of scenarios is very large or even infinite. Suppose further that we can generate a sample ${\displaystyle \xi ^{1},\xi ^{2},\dots ,\xi ^{N}}$ of ${\displaystyle N}$ replications of the random vector ${\displaystyle \xi }$. Usually the sample is assumed to be independent and identically distributed (i.i.d sample). Given a sample, the expectation function ${\displaystyle q(x)=E[Q(x,\xi )]}$ is approximated by the sample average

${\displaystyle {\hat {q}}_{N}(x)={\frac {1}{N}}\sum _{j=1}^{N}Q(x,\xi ^{j})}$

and consequently the first-stage problem is given by

${\displaystyle {\begin{array}{rlrrr}{\hat {g}}_{N}(x)=&\min \limits _{x\in \mathbb {R} ^{n}}&c^{T}x+{\frac {1}{N}}\sum _{j=1}^{N}Q(x,\xi ^{j})&\\&{\text{subject to}}&Ax&=&b\\&&x&\geq &0\end{array}}}$

This formulation is known as the Sample Average Approximation method. The SAA problem is a function of the considered sample and in that sense is random. For a given sample ${\displaystyle \xi ^{1},\xi ^{2},\dots ,\xi ^{N}}$ the SAA problem is of the same form as a two-stage stochastic linear programming problem with the scenarios ${\displaystyle \xi ^{j}}$., ${\displaystyle j=1,\dots ,N}$, each taken with the same probability ${\displaystyle p_{j}={\frac {1}{N}}}$.

Suppose the feasible set ${\displaystyle X}$ of the SAA problem is fixed, i.e., it is independent of the sample. Let ${\displaystyle \vartheta ^{*}}$ and ${\displaystyle S^{*}}$ be the optimal value and the set of optimal solutions, respectively, of the true problem and let ${\displaystyle {\hat {\vartheta }}_{N}}$ and ${\displaystyle {\hat {S}}_{N}}$ be the optimal value and the set of optimal solutions, respectively, of the SAA problem.

1. Let ${\displaystyle g:X\rightarrow \mathbb {R} }$ and ${\displaystyle {\hat {g}}_{N}:X\rightarrow \mathbb {R} }$ be a sequence of (deterministic) real valued functions. The following two properties are equivalent:
   • for any ${\displaystyle {\overline {x}}\in X}$ and any sequence ${\displaystyle \{x_{N}\}\subset X}$ converging to ${\displaystyle {\overline {x}}}$ it follows that ${\displaystyle {\hat {g}}_{N}(x_{N})}$ converges to ${\displaystyle g({\overline {x}})}$
   • the function ${\displaystyle g(\cdot )}$ is continuous on ${\displaystyle X}$ and ${\displaystyle {\hat {g}}_{N}(\cdot )}$ converges to ${\displaystyle g(\cdot )}$ uniformly on any compact subset of ${\displaystyle X}$
2. If the objective of the SAA problem ${\displaystyle {\hat {g}}_{N}(x)}$ converges to the true problem's objective ${\displaystyle g(x)}$ with probability 1, as ${\displaystyle N\rightarrow \infty }$, uniformly on the feasible set ${\displaystyle X}$. Then ${\displaystyle {\hat {\vartheta }}_{N}}$ converges to ${\displaystyle \vartheta ^{*}}$ with probability 1 as ${\displaystyle N\rightarrow \infty }$.
3. Suppose that there exists a compact set ${\displaystyle C\subset \mathbb {R} ^{n}}$ such that
   • the set ${\displaystyle S}$ of optimal solutions of the true problem is nonempty and is contained in ${\displaystyle C}$
   • the function ${\displaystyle g(x)}$ is finite valued and continuous on ${\displaystyle C}$
   • the sequence of functions ${\displaystyle {\hat {g}}_{N}(x)}$ converges to ${\displaystyle g(x)}$ with probability 1, as ${\displaystyle N\rightarrow \infty }$, uniformly in ${\displaystyle x\in C}$
   • for ${\displaystyle N}$ large enough the set ${\displaystyle {\hat {S}}_{N}}$ is nonempty and ${\displaystyle {\hat {S}}_{N}\subset C}$ with probability 1

then ${\displaystyle {\hat {\vartheta }}_{N}\rightarrow \vartheta ^{*}}$ and ${\displaystyle \mathbb {D} (S^{*},{\hat {S}}_{N})\rightarrow 0}$ with probability 1 as ${\displaystyle N\rightarrow \infty }$. Note that ${\displaystyle \mathbb {D} (A,B)}$ denotes the deviation of set ${\displaystyle A}$ from set ${\displaystyle B}$, defined as

${\displaystyle \mathbb {D} (A,B):=\sup _{x\in A}\{\inf _{x'\in B}\|x-x'\|\}}$

In some situations the feasible set ${\displaystyle X}$ of the SAA problem is estimated, then the corresponding SAA problem takes the form

${\displaystyle \min _{x\in X_{N}}{\hat {g}}_{N}(x)}$

where ${\displaystyle X_{N}}$ is a subset of ${\displaystyle \mathbb {R} ^{n}}$ depending on the sample and therefore is random. Nevertheless, consistency results for SAA estimators can still be derived under some additional assumptions:

1. Suppose that there exists a compact set ${\displaystyle C\subset \mathbb {R} ^{n}}$ such that
   • the set ${\displaystyle S}$ of optimal solutions of the true problem is nonempty and is contained in ${\displaystyle C}$
   • the function ${\displaystyle g(x)}$ is finite valued and continuous on ${\displaystyle C}$
   • the sequence of functions ${\displaystyle {\hat {g}}_{N}(x)}$ converges to ${\displaystyle g(x)}$ with probability 1, as ${\displaystyle N\rightarrow \infty }$, uniformly in ${\displaystyle x\in C}$
   • for ${\displaystyle N}$ large enough the set ${\displaystyle {\hat {S}}_{N}}$ is nonempty and ${\displaystyle {\hat {S}}_{N}\subset C}$ with probability 1
   • if ${\displaystyle x_{N}\in X_{N}}$ and ${\displaystyle x_{N}}$ converges with probability 1 to a point ${\displaystyle x}$, then ${\displaystyle x\in X}$
   • for some point ${\displaystyle x\in S^{*}}$ there exists a sequence ${\displaystyle x_{N}\in X_{N}}$ such that ${\displaystyle x_{N}\rightarrow x}$ with probability 1.

then ${\displaystyle {\hat {\vartheta }}_{N}\rightarrow \vartheta ^{*}}$ and ${\displaystyle \mathbb {D} (S^{*},{\hat {S}}_{N})\rightarrow 0}$ with probability 1 as ${\displaystyle N\rightarrow \infty }$.

Suppose the sample ${\displaystyle \xi ^{1},\dots ,\xi ^{N}}$ is i.i.d. and fix a point ${\displaystyle x\in X}$. Then the sample average estimator ${\displaystyle {\hat {g}}_{N}(x)}$, of ${\displaystyle g(x)}$, is unbiased and has variance ${\displaystyle {\frac {1}{N}}\sigma ^{2}(x)}$, where ${\displaystyle \sigma ^{2}(x):=Var[Q(x,\xi )]}$ is supposed to be finite. Moreover, by the central limit theorem we have that

${\displaystyle {\sqrt {N}}[{\hat {g}}_{N}-g(x)]{\xrightarrow {\mathcal {D}}}Y_{x}}$

where ${\displaystyle {\xrightarrow {\mathcal {D}}}}$ denotes convergence in distribution and ${\displaystyle Y_{x}}$ has a normal distribution with mean ${\displaystyle 0}$ and variance ${\displaystyle \sigma ^{2}(x)}$, written as ${\displaystyle {\mathcal {N}}(0,\sigma ^{2}(x))}$.

In other words, ${\displaystyle {\hat {g}}_{N}(x)}$ has asymptotically normal distribution, i.e., for large ${\displaystyle N}$, ${\displaystyle {\hat {g}}_{N}(x)}$ has approximately normal distribution with mean ${\displaystyle g(x)}$ and variance ${\displaystyle {\frac {1}{N}}\sigma ^{2}(x)}$. This leads to the following (approximate) ${\displaystyle 100(1-\alpha )}$% confidence interval for ${\displaystyle f(x)}$:

${\displaystyle \left[{\hat {g}}_{N}(x)-z_{\alpha /2}{\frac {{\hat {\sigma }}(x)}{\sqrt {N}}},{\hat {g}}_{N}(x)+z_{\alpha /2}{\frac {{\hat {\sigma }}(x)}{\sqrt {N}}}\right]}$

where ${\displaystyle z_{\alpha /2}:=\Phi ^{-1}(1-\alpha /2)}$ (here ${\displaystyle \Phi (\cdot )}$ denotes the cdf of the standard normal distribution) and

${\displaystyle {\hat {\sigma }}^{2}(x):={\frac {1}{N-1}}\sum _{j=1}^{N}\left[Q(x,\xi ^{j})-{\frac {1}{N}}\sum _{j=1}^{N}Q(x,\xi ^{j})\right]^{2}}$

is the sample variance estimate of ${\displaystyle \sigma ^{2}(x)}$. That is, the error of estimation of ${\displaystyle g(x)}$ is (stochastically) of order ${\displaystyle O({\sqrt {N}})}$.

Stochastic dynamic programming is frequently used to model animal behaviour in such fields as behavioural ecology.[8][9] Empirical tests of models of optimal foraging, life-history transitions such as fledging in birds and egg laying in parasitoid wasps have shown the value of this modelling technique in explaining the evolution of behavioural decision making. These models are typically many-staged, rather than two-staged.

Stochastic dynamic programming is a useful tool in understanding decision making under uncertainty. The accumulation of capital stock under uncertainty is one example; often it is used by resource economists to analyze bioeconomic problems[10] where the uncertainty enters in such as weather, etc.

The following is an example from finance of multi-stage stochastic programming. Suppose that at time ${\displaystyle t=0}$ we have initial capital ${\displaystyle W_{0}}$ to invest in ${\displaystyle n}$ assets. Suppose further that we are allowed to rebalance our portfolio at times ${\displaystyle t=1,\dots ,T-1}$ but without injecting additional cash into it. At each period ${\displaystyle t}$ we make a decision about redistributing the current wealth ${\displaystyle W_{t}}$ among the ${\displaystyle n}$ assets. Let ${\displaystyle x_{0}=(x_{10},\dots ,x_{n0})}$ be the initial amounts invested in the n assets. We require that each ${\displaystyle x_{i0}}$ is nonnegative and that the balance equation ${\displaystyle \sum _{i=1}^{n}x_{i0}=W_{0}}$ should hold.

Consider the total returns ${\displaystyle \xi _{t}=(\xi _{1t},\dots ,\xi _{nt})}$ for each period ${\displaystyle t=1,\dots ,T}$. This forms a vector-valued random process ${\displaystyle \xi _{1},\dots ,\xi _{T}}$. At time period ${\displaystyle t=1}$, we can rebalance the portfolio by specifying the amounts ${\displaystyle x_{1}=(x_{11},\dots ,x_{n1})}$ invested in the respective assets. At that time the returns in the first period have been realized so it is reasonable to use this information in the rebalancing decision. Thus, the second-stage decisions, at time ${\displaystyle t=1}$, are actually functions of realization of the random vector ${\displaystyle \xi _{1}}$, i.e., ${\displaystyle x_{1}=x_{1}(\xi _{1})}$. Similarly, at time ${\displaystyle t}$ the decision ${\displaystyle x_{t}=(x_{1t},\dots ,x_{nt})}$ is a function ${\displaystyle x_{t}=x_{t}(\xi _{[t]})}$ of the available information given by ${\displaystyle \xi _{[t]}=(\xi _{1},\dots ,\xi _{t})}$ the history of the random process up to time ${\displaystyle t}$. A sequence of functions ${\displaystyle x_{t}=x_{t}(\xi _{[t]})}$, ${\displaystyle t=0,\dots ,T-1}$, with ${\displaystyle x_{0}}$ being constant, defines an implementable policy of the decision process. It is said that such a policy is feasible if it satisfies the model constraints with probability 1, i.e., the nonnegativity constraints ${\displaystyle x_{it}(\xi _{[t]})\geq 0}$, ${\displaystyle i=1,\dots ,n}$, ${\displaystyle t=0,\dots ,T-1}$, and the balance of wealth constraints,

${\displaystyle \sum _{i=1}^{n}x_{it}(\xi _{[t]})=W_{t},}$

where in period ${\displaystyle t=1,\dots ,T}$ the wealth ${\displaystyle W_{t}}$ is given by

${\displaystyle W_{t}=\sum _{i=1}^{n}\xi _{it}x_{i,t-1}(\xi _{[t-1]}),}$

which depends on the realization of the random process and the decisions up to time ${\displaystyle t}$.

Suppose the objective is to maximize the expected utility of this wealth at the last period, that is, to consider the problem

${\displaystyle \max E[U(W_{T})].}$

This is a multistage stochastic programming problem, where stages are numbered from ${\displaystyle t=0}$ to ${\displaystyle t=T-1}$. Optimization is performed over all implementable and feasible policies. To complete the problem description one also needs to define the probability distribution of the random process ${\displaystyle \xi _{1},\dots ,\xi _{T}}$. This can be done in various ways. For example, one can construct a particular scenario tree defining time evolution of the process. If at every stage the random return of each asset is allowed to have two continuations, independent of other assets, then the total number of scenarios is ${\displaystyle 2^{nT}.}$

In order to write dynamic programming equations, consider the above multistage problem backward in time. At the last stage ${\displaystyle t=T-1}$, a realization ${\displaystyle \xi _{[T-1]}=(\xi _{1},\dots ,\xi _{T-1})}$ of the random process is known and ${\displaystyle x_{T-2}}$ has been chosen. Therefore, one needs to solve the following problem

${\displaystyle {\begin{array}{lrclr}\max \limits _{x_{T-1}}&E[U(W_{T})|\xi _{[T-1]}]&\\{\text{subject to}}&W_{T}&=&\sum _{i=1}^{n}\xi _{iT}x_{i,T-1}\\&\sum _{i=1}^{n}x_{i,T-1}&=&W_{T-1}\\&x_{T-1}&\geq &0\end{array}}}$

where ${\displaystyle E[U(W_{T})|\xi _{[T-1]}]}$ denotes the conditional expectation of ${\displaystyle U(W_{T})}$ given ${\displaystyle \xi _{[T-1]}}$. The optimal value of the above problem depends on ${\displaystyle W_{T-1}}$ and ${\displaystyle \xi _{[T-1]}}$ and is denoted ${\displaystyle Q_{T-1}(W_{T-1},\xi _{[T-1]})}$.

Similarly, at stages ${\displaystyle t=T-2,\dots ,1}$, one should solve the problem

${\displaystyle {\begin{array}{lrclr}\max \limits _{x_{t}}&E[Q_{t+1}(W_{t+1},\xi _{[t+1]})|\xi _{[t]}]&\\{\text{subject to}}&W_{t+1}&=&\sum _{i=1}^{n}\xi _{i,t+1}x_{i,t}\\&\sum _{i=1}^{n}x_{i,t}&=&W_{t}\\&x_{t}&\geq &0\end{array}}}$

whose optimal value is denoted by ${\displaystyle Q_{t}(W_{t},\xi _{[t]})}$. Finally, at stage ${\displaystyle t=0}$, one solves the problem

${\displaystyle {\begin{array}{lrclr}\max \limits _{x_{0}}&E[Q_{1}(W_{1},\xi _{[1]})]&\\{\text{subject to}}&W_{1}&=&\sum _{i=1}^{n}\xi _{i,1}x_{i0}\\&\sum _{i=1}^{n}x_{i0}&=&W_{0}\\&x_{0}&\geq &0\end{array}}}$

For a general distribution of the process ${\displaystyle \xi _{t}}$, it may be hard to solve these dynamic programming equations. The situation simplifies dramatically if the process ${\displaystyle \xi _{t}}$ is stagewise independent, i.e., ${\displaystyle \xi _{t}}$ is (stochastically) independent of ${\displaystyle \xi _{1},\dots ,\xi _{t-1}}$ for ${\displaystyle t=2,\dots ,T}$. In this case, the corresponding conditional expectations become unconditional expectations, and the function ${\displaystyle Q_{t}(W_{t})}$, ${\displaystyle t=1,\dots ,T-1}$ does not depend on ${\displaystyle \xi _{[t]}}$. That is, ${\displaystyle Q_{T-1}(W_{T-1})}$ is the optimal value of the problem

${\displaystyle {\begin{array}{lrclr}\max \limits _{x_{T-1}}&E[U(W_{T})]&\\{\text{subject to}}&W_{T}&=&\sum _{i=1}^{n}\xi _{iT}x_{i,T-1}\\&\sum _{i=1}^{n}x_{i,T-1}&=&W_{T-1}\\&x_{T-1}&\geq &0\end{array}}}$

and ${\displaystyle Q_{t}(W_{t})}$ is the optimal value of

${\displaystyle {\begin{array}{lrclr}\max \limits _{x_{t}}&E[Q_{t+1}(W_{t+1})]&\\{\text{subject to}}&W_{t+1}&=&\sum _{i=1}^{n}\xi _{i,t+1}x_{i,t}\\&\sum _{i=1}^{n}x_{i,t}&=&W_{t}\\&x_{t}&\geq &0\end{array}}}$

for ${\displaystyle t=T-2,\dots ,1}$.

All discrete stochastic programming problems can be represented with any algebraic modeling language, manually implementing explicit or implicit non-anticipativity to make sure the resulting model respects the structure of the information made available at each stage. An instance of an SP problem generated by a general modelling language tends to grow quite large (linearly in the number of scenarios), and its matrix loses the structure that is intrinsic to this class of problems, which could otherwise be exploited at solution time by specific decomposition algorithms. Extensions to modelling languages specifically designed for SP are starting to appear, see:

• AIMMS - supports the definition of SP problems
• EMP SP (Extended Mathematical Programming for Stochastic Programming) - a module of GAMS created to facilitate stochastic programming (includes keywords for parametric distributions, chance constraints and risk measures such as Value at risk and Expected shortfall).
• SAMPL - a set of extensions to AMPL specifically designed to express stochastic programs (includes syntax for chance constraints, integrated chance constraints and Robust Optimization problems)

They both can generate SMPS instance level format, which conveys in a non-redundant form the structure of the problem to the solver.