Quick start

Installation

StochasticPrograms is not yet registered and is therefore installed as follows:

pkg> add https://github.com/martinbiel/StochasticPrograms.jl

Afterwards, the functionality can be made available in a module or REPL through:

using StochasticPrograms

A simple stochastic program

To showcase the use of StochasticPrograms we will walk through a simple example. Consider the following stochastic program: (taken from Introduction to Stochastic Programming).

where

and the stochastic variable

takes on the value

with probability $0.4$ and

with probability $0.6$. In the following, we consider how to model, analyze, and solve this stochastic program using StochasticPrograms.

Scenario definition

First, we introduce a scenario type that can encompass the scenarios $\xi_1$ and $\xi_2$ above. This can be achieved conviently through the @scenario macro:

@scenario Simple = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end

Now, $\xi_1$ and $\xi_2$ can be created through:

ξ₁ = SimpleScenario(-24.0, -28.0, 500.0, 100.0, probability = 0.4)

SimpleScenario with probability 0.4

and

ξ₂ = SimpleScenario(-28.0, -32.0, 300.0, 300.0, probability = 0.6)

SimpleScenario with probability 0.6

Some useful functionality is automatically made available when introducing scenarios in this way. For example, we can check the discrete probability of a given scenario occuring:

probability(ξ₁)

0.4

Moreover, we can form the expected scenario out of a given set:

ξ̄ = expected([ξ₁, ξ₂])

Expected scenario of type SimpleScenario

Stochastic program definition

We are now ready to create a stochastic program based on the introduced scenario type. Optionally, we can also supply a capable MathProgBase solver that can be used internally when necessary. Consider:

using GLPKMathProgInterface

sp = StochasticProgram([ξ₁, ξ₂], solver = GLPKSolverLP())

Stochastic program with:
 * 2 scenarios of type SimpleScenario
 * undefined first stage
 * undefined second stage
Solver is GLPKInterfaceLP

The above command creates a stochastic program and preloads the two defined scenarios. The provided solver will be used internally when necessary. For clarity, we will still explicitly supply a solver when it is required. Now, we provide model recipes for the first and second stage of the example problem. The first stage is straightforward, and is defined using JuMP syntax inside a @first_stage block:

@first_stage sp = begin
    @variable(model, x₁ >= 40)
    @variable(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁ + x₂ <= 120)
end

Stochastic program with:
 * 2 scenarios of type SimpleScenario
 * 2 decision variables
 * undefined second stage
Solver is GLPKInterfaceLP

The recipe was immediately used to generate an instance of the first stage model. Next, we give a second stage recipe inside a @second_stage block:

@second_stage sp = begin
    @decision x₁ x₂
    ξ = scenario
    @variable(model, 0 <= y₁ <= ξ.d₁)
    @variable(model, 0 <= y₂ <= ξ.d₂)
    @objective(model, Min, ξ.q₁*y₁ + ξ.q₂*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end

Stochastic program with:
 * 2 scenarios of type SimpleScenario
 * 2 decision variables
 * 2 recourse variables
Solver is GLPKInterfaceLP

Every first stage variable that occurs in the second stage model is annotated with @decision at the beginning of the definition. Moreover, the scenario data is referenced through scenario. Instances of the defined scenario SimpleScenario will be injected to create instances of the second stage model. The second stage recipe is immediately used to generate second stage models for each preloaded scenario. Hence, the stochastic program definition is complete. We can now print the program and confirm that it indeed models the example recourse problem given above:

print(sp)

First-stage
==============
Min 100 x₁ + 150 x₂
Subject to
 x₁ + x₂ ≤ 120
 x₁ ≥ 40
 x₂ ≥ 20

Second-stage
==============
Subproblem 1:
Min -24 y₁ - 28 y₂
Subject to
 -60 x₁ + 6 y₁ + 10 y₂ ≤ 0
 -80 x₂ + 8 y₁ + 5 y₂ ≤ 0
 0 ≤ y₁ ≤ 500
 0 ≤ y₂ ≤ 100

Subproblem 2:
Min -28 y₁ - 32 y₂
Subject to
 6 y₁ + 10 y₂ - 60 x₁ ≤ 0
 8 y₁ + 5 y₂ - 80 x₂ ≤ 0
 0 ≤ y₁ ≤ 300
 0 ≤ y₂ ≤ 300

Deterministically equivalent problem

Since the example problem is small it is straightforward to work out the extended form:

which is also commonly referred to as the deterministically equivalent problem. This construct is available in StochasticPrograms through:

dep = DEP(sp)
print(dep)

Min 100 x₁ + 150 x₂ - 9.600000000000001 y₁_1 - 11.200000000000001 y₂_1 - 16.8 y₁_2 - 19.2 y₂_2
Subject to
 x₁ + x₂ ≤ 120
 6 y₁_1 + 10 y₂_1 - 60 x₁ ≤ 0
 8 y₁_1 + 5 y₂_1 - 80 x₂ ≤ 0
 6 y₁_2 + 10 y₂_2 - 60 x₁ ≤ 0
 8 y₁_2 + 5 y₂_2 - 80 x₂ ≤ 0
 x₁ ≥ 40
 x₂ ≥ 20
 0 ≤ y₁_1 ≤ 500
 0 ≤ y₂_1 ≤ 100
 0 ≤ y₁_2 ≤ 300
 0 ≤ y₂_2 ≤ 300

Evaluate decisions

With the stochastic program defined, we can now evaluate the performance of different first stage decisions. Consider the following first stage decision:

x = [40., 20.]

2-element Array{Float64,1}:
 40.0
 20.0

The expected result of taking this decision can be determined through:

evaluate_decision(sp, x, solver = GLPKSolverLP())

-470.39999999999964

The supplied solver is used to solve all available second stage models, with fixed first stage values. These outcome models can be built manually by supplying a scenario and the first stage decision.

print(outcome_model(sp, ξ₁, x))

Min -24 y₁ - 28 y₂
Subject to
 6 y₁ + 10 y₂ - 60 x₁ ≤ 0
 8 y₁ + 5 y₂ - 80 x₂ ≤ 0
 x₁ = 40
 x₂ = 20
 0 ≤ y₁ ≤ 500
 0 ≤ y₂ ≤ 100

Moreover, we can evaluate the result of the decision in a given scenario, i.e. solving a single outcome model, through:

evaluate_decision(sp, ξ₁, x, solver = GLPKSolverLP())

900.0

Optimal first stage decision

The optimal first stage decision is the decision that gives the best expected result over all available scenarios. This decision can be determined by solving the deterministically equivalent problem, by supplying a capable solver. Structure exploiting solvers are outlined in Structured solvers. In addition, it is possible to give a MathProgBase solver capable of solving linear programs. For example, we can solve sp with the GLPK solver as follows:

optimize!(sp, solver = GLPKSolverLP())

:Optimal

Internally, this generates and solves the extended form of sp. We can now inspect the optimal first stage decision through:

x_opt = optimal_decision(sp)

2-element Array{Float64,1}:
 46.66666666666667
 36.25

Moreover, the optimal value, i.e. the expected outcome of using the optimal decision, is acquired through:

optimal_value(sp)

-855.8333333333321

which of course coincides with the result of evaluating the optimal decision:

evaluate_decision(sp, x_opt, solver = GLPKSolverLP())

-855.8333333333321

This value is commonly referred to as the value of the recourse problem (VRP). We can also calculate it directly through:

VRP(sp, solver = GLPKSolverLP())

-855.8333333333321

Wait-and-see models

If we assume that we know what the actual outcome will be, we would be interested in the optimal course of action in that scenario. This is the concept of wait-and-see models. For example if $ξ₁$ is believed to be the actual outcome, we can define a wait-and-see model as follows:

ws = WS(sp, ξ₁)
print(ws)

Min 100 x₁ + 150 x₂ - 24 y₁ - 28 y₂
Subject to
 x₁ + x₂ ≤ 120
 6 y₁ + 10 y₂ - 60 x₁ ≤ 0
 8 y₁ + 5 y₂ - 80 x₂ ≤ 0
 x₁ ≥ 40
 x₂ ≥ 20
 0 ≤ y₁ ≤ 500
 0 ≤ y₂ ≤ 100

The optimal first stage decision in this scenario can be determined through:

x₁ = WS_decision(sp, ξ₁, solver = GLPKSolverLP())

2-element Array{Float64,1}:
 40.0
 29.583333333333336

We can evaluate this decision:

evaluate_decision(sp, x₁, solver = GLPKSolverLP())

-762.5

The outcome is of course worse than taking the optimal decision. However, it would perform better if $ξ₁$ is the actual outcome:

evaluate_decision(sp, ξ₁, x₁, solver = GLPKSolverLP())

37.5

as compared to:

evaluate_decision(sp, ξ₁, x_opt, solver = GLPKSolverLP())

104.16666666666788

Another important concept is the wait-and-see model corresponding to the expected future scenario. This is referred to as the expected value problem and can be generated through:

evp = EVP(sp)
print(evp)

Min 100 x₁ + 150 x₂ - 26.400000000000002 y₁ - 30.4 y₂
Subject to
 x₁ + x₂ ≤ 120
 6 y₁ + 10 y₂ - 60 x₁ ≤ 0
 8 y₁ + 5 y₂ - 80 x₂ ≤ 0
 x₁ ≥ 40
 x₂ ≥ 20
 0 ≤ y₁ ≤ 380
 0 ≤ y₂ ≤ 220

Internally, this generates the expected scenario out of the available scenarios and forms the respective wait-and-see model. The optimal first stage decision associated with the expected value problem is conviently determined using

x̄ = EVP_decision(sp, solver = GLPKSolverLP())

2-element Array{Float64,1}:
 71.45833333333334
 48.54166666666667

Again, we can evaluate this decision:

evaluate_decision(sp, x̄, solver = GLPKSolverLP())

-568.9166666666679

This value is often referred to as the expected result of using the expected value solution (EEV), and is also available through:

EEV(sp, solver = GLPKSolverLP())

-568.9166666666679

Stochastic performance

Finally, we consider some performance measures of the defined model. The expected value of perfect information is the difference between the value of the recourse problem and the expected result of having perfect knowledge. In other words, it involes solving the recourse problem as well as every wait-and-see model that can be formed from the available scenarios. We calculate it as follows:

EVPI(sp, solver = GLPKSolverLP())

662.9166666666679

The resulting value indicates the expected gain of having perfect information about future scenarios. Another concept is the value of the stochastic solution, which is the difference between the value of the recourse problem and the EEV. We calculate it as follows:

VSS(sp, solver = GLPKSolverLP())

286.91666666666424

The resulting value indicates the gain of including uncertainty in the model formulation.