# Solving Stochastic Models

Solving stochastic models

After generating stochastic event data and a scenario tree, you can generate and solve the stochastic model by using methods from the GMP library discussed in Implementing Advanced Algorithms for Mathematical Programs. AIMMS supports two methods for solving a stochastic model:

by solving its

*deterministic equivalent*, orfor purely linear mathematical programs only, through the

*stochastic Benders algorithm*, orusing the Benders decomposition algorithm in CPLEX 12.7 or higher.

Both Benders algorithms will decompose the stochastic model into
multiple smaller models, and thus is better suited to solve stochastic
models where the deterministic equivalent, either by the size of the
deterministic model or because of a huge number of scenarios, becomes
too big or time-consuming to solve at once. The Benders decomposition
algorithm in CPLEX can be used to solve stochastic models with integer
variables, as long as all integer variables are assigned to the first
stage. For more information see the CPLEX option `Benders_strategy`

.

## Generating and Solving the Deterministic Equivalent

Generating a stochastic model

The method for generating a stochastic model for a
`MathematicalProgram`

*MP* is

`GMP::Instance::GenerateStochasticProgram`

(*MP*,*StochasticParameters*,*StochasticVariables*,*Scenarios*,*ScenarioProbability*,*ScenarioTreeMap*,*DeterministicScenarioName*[,*GenerationMode*][,*Name*])

The function returns an element into the set
`AllGeneratedMathematicalPrograms`

. This generated math program
instance contains a memory-efficient representation of the technology
matrix of the stochastic model and the stochastic event data, and can be
used to create a deterministic equivalent of the stochastic model, as
well as the submodels necessary for a stochastic Benders approach.

Specifying stochastic identifiers

Through the arguments *StochasticParameters* and *StochasticVariables*
you indicate to AIMMS which stochastic parameters and variables you want
to take into consideration when generating this stochastic model. These
arguments must be subsets of the predefined sets
`AllStochasticParameters`

and `AllStochasticVariables`

,
respectively. You may want to use real subsets, for instance, when your
AIMMS project contains multiple stochastic models, each referring only
to a subset of the stochastic parameters and variables.

Specifying scenarios

Through the *Scenarios*, *ScenarioProbability* and *ScenarioTreeMap*
arguments you specify the set of scenarios, their probabilities and the
mapping defining the scenario tree for which you want to generate the
stochastic model to AIMMS. Through the string argument
*DeterministicScenarioName*, you supply the name of the artificial
element that AIMMS will add to the predefined set
`AllStochasticScenarios`

(if not created already), and use to store
the solution of non-stochastic variables in their respective
.Stochastic suffices as explained in Stochastic Parameters and Variables.

Enforcing non-anticipativity constraints

Using the *GenerationMode* argument you can specify whether you want
AIMMS to explicitly add the non-anticipativity constraints to your
stochastic model, or whether you want non-anticipativity to be enforced
implicitly by substituting the representative scenario for every
non-representative scenario at every stage. *GenerationMode* is an
element parameter into the predefined set
`AllStochasticGenerationModes`

, with possible values

`'CreateNonAnticipativityConstraints'`

, and`'SubstituteStochasticVariables'`

(the default value).

Name argument

With the optional *Name* argument you can explicitly specify a name for
the generated mathematical program. If you do not choose a name, AIMMS
will use the name of the underlying `MathematicalProgram`

as the name
of the generated mathematical program as well. Please note, that AIMMS
will also use this name as the default name for solving the
deterministic model. Therefore, if you do not want the generated
mathematical program of the deterministic model to be deleted, then you
have to choose a non-default name.

Solving the deterministic equivalent of a stochastic model

You can solve a stochastic model by using the regular GMP procedure

`GMP::Instance::Solve`

(*gmp*)

By applying this function to a generated mathematical program associated
with a stochastic model, AIMMS will create the deterministic equivalent
and pass it to the appropriate LP/MIP solver. The
`GMP::Instance::Solve`

method is discussed in full detail in
Managing Generated Mathematical Program Instances.

Changing the model input

Note that, when you adjust the scenario tree map, the stochastic data,
the scenario probabilities, or the value of the `Stage`

attribute of
some variables after you generated the stochastic model, you should
regenerate the stochastic model again to reflect these changes.

Example

Consider the following call to
`GMP::Instance::GenerateStochasticProgram`

```
GMP::Instance::GenerateStochasticProgram(
TransportModel, AllStochasticParameters, AllStochasticVariables,
MyScenarios, MyScenarioProbability, MyScenarioTreeMap,
"TransportModel", 'SubstituteStochasticVariables', "StochasticTransportModel");
```

After solving the generated stochastic model, its solution will be
stored as follows, where `sc`

is an index into `MyScenarios`

the per-scenario solution of a stochastic variable

`Transport(i,j)`

will be stored in`Transport.Stochastic(sc,i,j)`

,the deterministic solution of a non-stochastic variable

`InitialStock(i)`

will be stored in`InitialStock.Stochastic('TransportModel',i)`

,the weighted objective value for the objective variable

`TotalCost`

will be stored in`TotalObjective.Stochastic('TransportModel')`

, while the contribution by every scenario is available through`TotalCost.Stochastic(sc)`

.

## Using the Stochastic Benders Algorithm

Using the stochastic Benders algorithm

Instead of solving the deterministic equivalent of a stochastic model,
AIMMS also allows you to solve *linear* stochastic models using a
stochastic Benders algorithm. The stochastic Benders algorithm is based
on a reformulation of the original model as a sequence of models
outlined below. The solution of the original model can be achieved by
solving the sequence of models iteratively until a terminating condition
is reached. A more detailed discussion of the stochastic Benders
algorithm can be found in [DT98] or [Alt03].

Definitions

All nodes in the scenario tree are numbered starting at 1 (the root node).

Convention

In the algorithmic outline below we identify the problem names with their associated solutions. That is, if a problem is, for instance, identified as \(F_i(x_{p_i})\), we will also use this name to denote its solution in other sub- problems.

The original model

The nested Benders algorithm can be used for problems of the form

A reformulation as a sequence of models

This problem corresponds to the following sequence of problems. For node \(i = 1\), the problem \(F_1\) is defined as

For all other nodes \(i\in N_t\) in stage \(t\in T\backslash\{1\}\), the problem \(F_i(x_{p_i})\) is defined as follows (note that \(\sum_{j\in I_i} q_j = q_i\))

For the leaf nodes in the scenario tree, the term \(\sum_{j\in I_i}\frac{q_j}{q_i}F_j(x_i)\) is omitted.

Formulated differently

If we now introduce an upper bound \(\theta_i\) to replace the term \(\sum_{j\in I_i}\frac{q_j}{q_i}F_j(x_i)\), we can rewrite the subproblem \(F_i(x_{p_i})\) as

Because of the linear nature of the original problem, the terms \(\sum_{j\in I_i}\frac{q_j}{q_i}F_j(x_i)\) are piecewise linear and convex. Therefore there exists an (a priori unknown) set of equations

that describes such a term and for which

Moreover, we are only interested in those \(x_i\) such that \(F_j(x_i)\) are feasible for all \(j\in I_i\). This requirement can be enforced by an (a priori unknown) set of constraints

By substituting these constraints we obtain the following reformulation of problem \(F_i(x_{p_i})\)

The relaxed master problem

The actual problem that is solved at node \(i\) is the following relaxed master problem \(\tilde{N}_i(x_{p_i})\) defined as follows:

At the start of the Benders algorithm \(r_i\) and \(s_i\) will be 0 for all \(i\in N\). The constraints \(D^l_i x_i + \theta_i\geq d^l_i\) are optimality cuts obtained from the children. That is, if \(\tilde{N}_j(x_i)\) is feasible for all \(j\in I_i\) (but not optimal) then an optimality cut is added to \(\tilde{N}_i(x_{p_i})\). The optimality cut is constructed by using a combination of the dual solutions of \(\tilde{N}_j(x_i)\) for all \(j\in I_i\). Adding an optimality cut does not make a feasible relaxed master problem infeasible. The Benders algorithm fails if one of the subproblems is unbounded. This can be avoided by giving all variables, except the objective variable, finite bounds.

Adding feasibility cuts

The constraints \(E^l_i x_i\geq e^l_i\) are feasibility cuts obtained from a child. If some child problem \(\tilde{N}_j(x_i)\) is not feasible then the following problem \(\tilde{E}_j(x_i)\) is solved

This feasibility problem can only be formulated for linear problems, is always feasible, and bounded from below by 0. Its dual solution is used to construct a new feasibility constraint for \(\tilde{N}_i(x_{p_i})\). Note that node \(j\) in its turn obtains optimality and/or feasibility cuts from its children for \(\tilde{N}_j(x_i)\) and \(\tilde{E}_j(x_i)\), unless \(j\) refers to a leaf node.

Terminating condition

If \((x_i, \theta_i)\) is an optimal solution of \(\tilde{N}_i(x_{p_i})\) and

then \((x_i, \theta_i)\) is an optimal solution of \(F_i(x_{p_i})\). If this holds for all non-leaf nodes then we have found an optimal solution of our original problem. For the leaf nodes, \(x_i\) only needs to be an optimal solution of \(\tilde{N}_i(x_{p_i})\).

Implementation in AIMMS

The stochastic Benders algorithm outlined above is implemented in AIMMS as a system module that you can include into your model, together with a number of supporting functions in the GMP library to perform a number of algorithmic steps that cannot be performed in the AIMMS language itself, for instance, to actually generate the stochastic sub-problems, and to generate feasibility and optimality cuts.

Adding the module

You can add the system module implementing the stochastic Benders
algorithm to your model through the **Settings-Install System Module…**
menu. By selecting the **Stochastic Decomposition Module** in the
**Install System Module** dialog box, AIMMS will add this system module
to your model.

Using the stochastic Benders module

The main procedure for using the stochastic Benders algorithm is
`DoStochasticDecomposition`

. Its inputs are:

a stochastic GMP,

the set of stages to consider, and

the set of scenarios to consider.

The procedure implements the algorithm outlined above. The supporting GMP functions for the stochastic Benders algorithm are described in Supporting Functions for Stochastic Programs.

Modifying the algorithm

Because the stochastic Benders algorithm is written in the AIMMS
language, you have complete freedom to modify the algorithm in order to
tune it for your stochastic programs. Also, the basic algorithm solves
all sub-problems sequentially. If your AIMMS license permits parallel
solver sessions, you can also speed up the algorithm by solving multiple
sub-problems in parallel using the GMP function
`GMP::SolverSession::AsynchronousExecute`

.