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, or

  • for purely linear mathematical programs only, through the stochastic Benders algorithm, or

  • using 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

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

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).

\[\begin{split}\begin{align} & \textbf{Indices:} \\ &&& \text{$i$} & & \text{index for the set of nodes $N$} \\ &&& \text{$t$} & & \text{index for the set of stages $T$} \\ & \textbf{Parameters:} \\ &&& \text{$q_i$} & & \text{probability belonging to node $i$} \\ &&& \text{$p_i$} & & \text{parent of node $i$} \\ & \textbf{Sets:} \\ &&& \text{$I_i$} & & \text{set with children of node $i$} \\ &&& \text{$N_t$} & & \text{set of nodes belonging to stage $t$} \end{align}\end{split}\]

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

\[\begin{split}\begin{align} & \text{minimize} & & \sum_{t\in T}\sum_{i\in N_t} q_i c_i^T x_i \\ & \text{subject to} & & W_1 x_1 = h_1 & & \\ &&& A_i x_{p_i} + W_i x_i = h_i & & \forall i\in N_t, t\in T\backslash\{1\} \\ &&& x_i \geq 0 & & \forall i\in N_t, t\in T \\ \end{align}\end{split}\]

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

\[\begin{split}\begin{align} & \text{minimize} & & c_1^Tx_1 + \sum_{j\in I_1} q_i F_j(x_1) \\ & \text{subject to} & & W_1x_1 = h_1 & & \\ &&& x_1 \geq 0 & & \\ \end{align}\end{split}\]

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\))

\[\begin{split}\begin{align} & \text{minimize} & & c_i^Tx_i + \sum_{j\in I_i} \frac{q_j}{q_i} F_j(x_i) \\ & \text{subject to} & & W_ix_i = h_i - A_ix_{p_i} & & \\ &&& x_i \geq 0 & & \\ \end{align}\end{split}\]

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

\[\begin{split}\begin{align} & \text{minimize} & & c_i^Tx_i + \theta_i \\ & \text{subject to} & & W_ix_i = h_i - A_ix_{p_i} & & \\ &&& \theta_i \geq \sum_{j\in I_i}\frac{q_j}{q_i}F_j(x_i) & & \\ &&& x_i \geq 0 & & \\ \end{align}\end{split}\]

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

\[D^l_i x_i = d^l_i\]

that describes such a term and for which

\[D^l_i x_i + \theta_i \geq d^l_i.\]

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

\[E^l_i x_i \geq e^l_i.\]

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

\[\begin{split}\begin{align} & \text{minimize} & & c_i^Tx_i + \theta_i \\ & \text{subject to} & & W_ix_i = h_i - A_ix_{p_i} & & \\ &&& D^l_i x_i + \theta_i \geq d^l_i & & \forall l\in 1,\dots,R_i \\ &&& E^l_i x_i \geq e^l_i & & \forall l\in 1,\dots,S_i \\ &&& x_i \geq 0 & & \\ \end{align}\end{split}\]

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:

\[\begin{split}\begin{align} & \text{minimize} & & c_i^Tx_i + \theta_i \\ & \text{subject to} & & W_ix_i = h_i - A_ix_{p_i} & & \\ &&& D^l_i x_i + \theta_i \geq d^l_i & & \forall l\in 1,\dots,r_i \\ &&& E^l_i x_i \geq e^l_i & & \forall l\in 1,\dots,s_i \\ &&& x_i \geq 0 & & \\ \end{align}\end{split}\]

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

\[\begin{split}\begin{align} & \text{minimize} & & w_j = e^T u^+_j + e^T u^-_j \\ & \text{subject to} & & W_jx_j + I u^+_j - I u^-_j = h_j - A_jx_i & & \\ &&& E^l_j x_j \geq e^l_j & & \forall l\in 1,\dots,s_j \\ &&& x_j \geq 0 & & \\ &&& u^+_j \geq 0 & & \\ &&& u^-_j \geq 0 & & \\ \end{align}\end{split}\]

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

\[\theta_i \geq \tilde{N}_i(x_{p_i})\]

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.