Stochastic data

Decoupling data design and model design is a fundamental principle in StochasticPrograms. This decoupling is achieved through data injection. By data we mean parameters in an optimization problem. In StochasticPrograms, this data is either deterministic and related to a specific stage, or uncertain and related to a specific scenario.

Stage data

Stage data is related to parameters that always appear in the first or second stage of a stochastic program. These parameters are deterministic and are the same across all scenarios. Stage data must be supplied when creating a stochastic program, through specialized constructors. However, the data can later be mutated without having to construct a new stochastic program instance. Any stage related data can then be accessed through the reserved keyword stage in @first_stage and @second_stage blocks. To showcase, we consider a minimal stochastic program:

\[\DeclareMathOperator*{\maximize}{maximize} \begin{aligned} \maximize_{x \in \mathbb{R}} & \quad x + \operatorname{\mathbb{E}}_{\omega} \left[Q(x, \xi(\omega))\right] \\ \text{s.t.} & \quad l_1 \leq x \leq u_1 \end{aligned}\]

where

\[\begin{aligned} Q(x, \xi(\omega)) = \max_{y \in \mathbb{R}} & \quad q_{\omega} y \\ \text{s.t.} & \quad y + x \leq U \\ & \quad l_2 \leq y \leq u_2 \end{aligned}\]

and the stochastic variable

\[ \xi(\omega) = q_{\omega}\]

takes on the value $1$ or $-1$ with equal probability. Here, the first stage contains the two parameters: $l_1$ and $u_1$. The second stage contains the three scenario-independent parameters: $U$, $l_2$, and $u_2$. The following defines this problem in StochasticPrograms, with some chosen parameter values:

using StochasticPrograms
using GLPKMathProgInterface

@scenario Simple = begin
    q::Float64
end

ξ₁ = SimpleScenario(1., probability = 0.5)
ξ₂ = SimpleScenario(-1., probability = 0.5)

l₁ = -1.
u₁ = 1.

U = 2.
l₂ = -1.
u₂ = 1.

sp = StochasticProgram((l₁,u₁), (U,l₂,u₂), [ξ₁,ξ₂])

@first_stage sp = begin
    l₁, u₁ = stage
    @variable(model, l₁ <= x <= u₁)
    @objective(model, Max, x)
end

@second_stage sp = begin
    @decision x
    U, l₂, u₂ = stage
    ξ = scenario
    @variable(model, l₂ <= y <= u₂)
    @objective(model, Max, ξ.q*y)
    @constraint(model, y + x <= U)
end

print(sp)

print("VRP = $(VRP(sp, solver = GLPKSolverLP()))")

First-stage
==============
Max x
Subject to
 -1 ≤ x ≤ 1

Second-stage
==============
Subproblem 1:
Max y
Subject to
 x + y ≤ 2
 -1 ≤ y ≤ 1

Subproblem 2:
Max -y
Subject to
 y + x ≤ 2
 -1 ≤ y ≤ 1

VRP = 2.0

Now, we can investigate the impact of the stage parameters by changing them slightly and regenerate the problem:

l₁ = -2.
u₁ = 2.

U = 2.
l₂ = -0.5
u₂ = 0.5

set_first_stage_data!(sp, (l₁,u₁))
set_second_stage_data!(sp, (U,l₂,u₂))

generate!(sp) # Regenerate problem

print(sp)

print("VRP = $(VRP(sp, solver = GLPKSolverLP()))")

First-stage
==============
Max x
Subject to
 -2 ≤ x ≤ 2

Second-stage
==============
Subproblem 1:
Max y
Subject to
 y + x ≤ 2
 -0.5 ≤ y ≤ 0.5

Subproblem 2:
Max -y
Subject to
 y + x ≤ 2
 -0.5 ≤ y ≤ 0.5

VRP = 2.25

Scenario data

Any uncertain parameter in the second stage of a stochastic program should be included in some predefined AbstractScenario type. Hence, all uncertain parameters in a stochastic program must be identified before defining the models. In brief, StochasticPrograms demands two functions from this abstraction. The discrete probability of a given AbstractScenario occurring should be returned from probability. Also, the expected scenario out of a collection of given AbstractScenarios should be returned by expected. StochasticPrograms provides a convenience macro, @scenario, for creating scenario types that adhere to this abstraction. If the identified uncertain parameters is a collection of numerical values, this is the recommended way to define the required scenario type.

using StochasticPrograms

@scenario Example = begin
    X::Float64
end

s₁ = ExampleScenario(1., probability = 0.5)
s₂ = ExampleScenario(5., probability = 0.5)

println("Probability of s₁: $(probability(s₁))")

s = expected([s₁, s₂])

println("Expectation over s₁ and s₂: $s")
println("Expectated X: $(s.scenario.X)")

Probability of s₁: 0.5
Expectation over s₁ and s₂: Expected scenario of type ExampleScenario
Expectated X: 3.0

Here, all the required operations are correctly defined.

There are some caveats to note. First, the autogenerated requires an additive zero element of the introduced scenario type. For simple numeric types this is autogenerated as well. However, say that we want to extend the above scenario with some vector parameter of size 2:

using StochasticPrograms

@scenario Example = begin
    X::Float64
    Y::Vector{Float64}
end

┌ Warning: Zero not defined for Array{Float64,1}. Cannot generate zero function.
└ @ StochasticPrograms ~/build/martinbiel/StochasticPrograms.jl/src/methods/util.jl:178
┌ Warning: The scenario type Example was not defined. A user-provided implementation
│
│     function zero(::Type{{ExampleScenario})
│         ...
│     end
│
│ is required.
└ @ Main.##ex-#470 ~/build/martinbiel/StochasticPrograms.jl/src/methods/creation.jl:138

In this case, we must provide an implementation of zero using @zero:

using StochasticPrograms

@scenario Example = begin
    X::Float64
    Y::Vector{Float64}

    @zero begin
        return Example(0.0, [0.0, 0.0])
    end
end

s₁ = ExampleScenario(1., ones(2), probability = 0.5)
s₂ = ExampleScenario(5., -ones(2), probability = 0.5)

println("Probability of s₁: $(probability(s₁))")

s = expected([s₁, s₂])

println("Expectation over s₁ and s₂: $s")
println("Expectated X: $(s.scenario.X)")
println("Expectated Y: $(s.scenario.Y)")

Probability of s₁: 0.5
Expectation over s₁ and s₂: Expected scenario of type ExampleScenario
Expectated X: 3.0
Expectated Y: [0.0, 0.0]

Another caveat is that the expected function can only be auto generated for fields that support addition and scalar multiplication with Float64. Consider:

using StochasticPrograms

@scenario Example = begin
    X::Float64
    Y::Vector{Float64}
    Z::Int

    @zero begin
        return Example(0.0, [0.0, 0.0], 0)
    end
end

┌ Warning: Scalar multiplication with Float64 not defined for Int64. Cannot generate expectation function.
└ @ StochasticPrograms ~/build/martinbiel/StochasticPrograms.jl/src/methods/util.jl:192
┌ Warning: The scenario type Example was not defined. A user-provided implementation
│
│     function expected(scenarios::Vector{ExampleScenario})
│         ...
│     end
│
│ is required.
└ @ Main.##ex-#472 ~/build/martinbiel/StochasticPrograms.jl/src/methods/creation.jl:141

Again, the solution is to provide an implementation of expected, this time using @expectation:

using StochasticPrograms

@scenario Example = begin
    X::Float64
    Y::Vector{Float64}
    Z::Int

    @zero begin
        return Example(0.0, [0.0, 0.0], 0)
    end

    @expectation begin
        X = sum([probability(s)*s.X for s in scenarios])
        Y = sum([probability(s)*s.Y for s in scenarios])
        Z = sum([round(Int, probability(s)*s.Z) for s in scenarios])
        return Example(X, Y, Z)
    end
end

s₁ = ExampleScenario(1., ones(2), 1, probability = 0.5)
s₂ = ExampleScenario(5., -ones(2), -1, probability = 0.5)

println("Probability of s₁: $(probability(s₁))")

s = expected([s₁, s₂])

println("Expectation over s₁ and s₂: $s")
println("Expectated X: $(s.scenario.X)")
println("Expectated Y: $(s.scenario.Y)")
println("Expectated Z: $(s.scenario.Z)")

Probability of s₁: 0.5
Expectation over s₁ and s₂: Expected scenario of type ExampleScenario
Expectated X: 3.0
Expectated Y: [0.0, 0.0]
Expectated Z: 0

For most problems, @scenario will probably be adequate. Otherwise consider defining Custom scenarios.

Sampling

Typically, we do not have exact knowledge of all possible future scenarios. However, we often have access to some model of the uncertainty. For example, scenarios could originate from:

A stochastic variable with known distribution
A time series fitted to data
A nerual network prediction

Even if the exact scenario distribution is unknown, or not all possible scenarios are available, we can still formulate a stochastic program that approximates the model we wish to formulate. This is achieved through a technique called sampled average approximation, which is based on sampling. The idea is to sample a large number $n$ of scenarios with equal probability $\frac{1}{n}$ and then use them to generate and solve a stochastic program. By the law of large numbers, the result will converge with probability $1$ to the "true" solution with increasing $n$.

StochasticPrograms accepts AbstractSampler objects in place of AbstractScenario. However, an AbstractSampler is always linked to some underlying AbstractScenario type, which is reflected in the resulting stochastic program as well. Samplers are conviniently created using @sampler. We can define a simple scenario type and a simple sampler as follows:

using StochasticPrograms

@scenario Example = begin
    ξ::Float64
end

@sampler Example = begin
    w::Float64

    Example(w::AbstractFloat) = new(w)

    @sample begin
        w = sampler.w
        return ExampleScenario(w*randn(), probability = rand())
    end
end

This creates a new AbstractSampler type called ExampleSampler, which samples ExampleScenarios. Now, we can create a sampler object and sample a scenario

sampler = ExampleSampler(2.)

s = sampler()

println(s)
println("ξ: $(s.ξ)")

ExampleScenario with probability 0.34651701419196046
ξ: 0.5945759690709232

It is possible to create other sampler objects for the ExampleScenario, by providing a new unique name:

@sampler Another Example = begin
    w::Float64
    d::Float64

    Another(w::AbstractFloat, d::AbstractFloat) = new(w, d)

    @sample begin
        w = sampler.w
        d = sampler.d
        return ExampleScenario(w*randn() + d, probability = rand())
    end
end

another = AnotherSampler(2., 6.)

s = another()

println(s)
println("ξ: $(s.ξ)")

ExampleScenario with probability 0.00790928339056074
ξ: 4.804731046543537

Now, lets use the first sampler to create a stochastic program:

sp = StochasticProgram(sampler)

@first_stage sp = begin
    @variable(model, x >= 0)
    @objective(model, Min, x)
end

@second_stage sp = begin
    @decision x
    ξ = scenario.ξ
    @variable(model, y)
    @objective(model, Min, y)
    @constraint(model, y + x == ξ)
end

Stochastic program with:
 * 0 scenarios of type ExampleScenario
 * 1 decision variable
 * 0 recourse variables
Solver is default solver

Now, we can sample $5$ scenarios to generate $5$ subproblems:

sample!(sp, 5)

Stochastic program with:
 * 5 scenarios of type ExampleScenario
 * 1 decision variable
 * 1 recourse variable
Solver is default solver

Printing yields:

print(sp)

First-stage
==============
Min x
Subject to
 x ≥ 0

Second-stage
==============
Subproblem 1:
Min y
Subject to
 y + x = -1.678053708777528
 y

Subproblem 2:
Min y
Subject to
 x + y = -4.534172697601061
 y

Subproblem 3:
Min y
Subject to
 x + y = 1.1674165751375571
 y

Subproblem 4:
Min y
Subject to
 x + y = -1.0446735148430168
 y

Subproblem 5:
Min y
Subject to
 x + y = -1.3873072876077712
 y

Sampled stochastic programs are solved as usual:

using GLPKMathProgInterface

optimize!(sp, solver = GLPKSolverLP())

println("optimal decision: $(optimal_decision(sp))")
println("optimal value: $(optimal_value(sp))")

optimal decision: [0.0]
optimal value: -1.4953581267383638

SSA is a shorthand for the above sequence of commands, which also accepts another sampler object over the same scenario type. For example:

using GLPKMathProgInterface

res = SSA(sp, another, 5, solver = GLPKSolverLP())

println("optimal decision: $(optimal_decision(sp))")
println("optimal value: $res")

optimal decision: [5.86646]
optimal value: 4.919457996396745

The quality of the model can be checked in different ways. One indicator is:

VSS(sp, solver = GLPKSolverLP())

-8.881784197001252e-16

Another is acquired by evaluating the optimal decision on a larger number of sampled scenarios:

x = optimal_decision(sp)
sample!(sp, 10000)
evaluate_decision(sp, x, solver = GLPKSolverLP())

0.00972763180220948

Again, if the functionality offered by @sampler is not adequate, consider Custom scenarios.

Custom scenarios

More complex scenario designs are probably not implementable using @scenario. However, it is still possible to create a custom scenario type as long as:

The type is a subtype of AbstractScenario
The type implements probability
The type implements expected, which should return an additive zero element if given an empty array

The restriction on expected is there to support taking expectations in a distributed environment.

As an example, consider the following generalized stochastic program:

\[\DeclareMathOperator*{\minimize}{minimize} \begin{aligned} \minimize_{x \in \mathbb{R}} & \quad \operatorname{\mathbb{E}}_{\omega} \left[(x - \xi(\omega))^2\right] \\ \end{aligned}\]

where $\xi(\omega)$ is exponentially distributed. We will skip the mathematical details here and just take for granted that the optimizer to the above problem is the mean of the exponential distribution. We will try to approximately solve this problem using sample average approximation. First, lets try to introduce a custom discrete scenario type that models a stochastic variable with a continuous probability distribution. Consider the following implementation:

using StochasticPrograms
using Distributions

struct DistributionScenario{D <: UnivariateDistribution} <: AbstractScenario
    probability::Probability
    distribution::D
    ξ::Float64

    function DistributionScenario(distribution::UnivariateDistribution, val::AbstractFloat)
        return new{typeof(distribution)}(Probability(pdf(distribution, val)), distribution, Float64(val))
    end
end

function StochasticPrograms.expected(scenarios::Vector{<:DistributionScenario{D}}) where D <: UnivariateDistribution
    isempty(scenarios) && return DistributionScenario(D(), 0.0)
    distribution = scenarios[1].distribution
    return ExpectedScenario(DistributionScenario(distribution, mean(distribution)))
end

The fallback probability method is viable as long as the scenario type contains a Probability field named probability. The implementation of expected is somewhat unconventional as it returns the mean of the distribution regardless of how many scenarios are given.

We are also free to define custom sampler objects, as long as:

The sampler type is a subtype of AbstractSampler
The sampler type implements a functor call that performs the sampling

We can implement a sampler that generates exponentially distributed scenarios as follows:

struct ExponentialSampler <: AbstractSampler{DistributionScenario{Exponential{Float64}}}
    distribution::Exponential

    ExponentialSampler(θ::AbstractFloat) = new(Exponential(θ))
end

function (sampler::ExponentialSampler)()
    ξ = rand(sampler.distribution)
    return DistributionScenario(sampler.distribution, ξ)
end

Now, lets attempt to define the generalized stochastic program using the available modeling tools:

using Ipopt

sampler = ExponentialSampler(2.)
sp = StochasticProgram(sampler)

@first_stage sp = begin
    @variable(model, x)
end

@second_stage sp = begin
    @decision x
    ξ = scenario.ξ
    @variable(model, y)
    @constraint(model, y == (x - ξ)^2)
    @objective(model, Min, y)
end

Stochastic program with:
 * 0 scenarios of type DistributionScenario
 * 1 decision variable
 * 0 recourse variables
Solver is default solver

The mean of the given exponential distribution is $2.0$, which is the optimal solution to the general problem. Now, lets sample 1000 exponentially distributed numbers with equal probability:

StochasticPrograms.sample!(sp, 1000) # Sample 1000 exponentially distributed scenarios (qualified call due to name clash with Distributions.jl)

Stochastic program with:
 * 1000 scenarios of type DistributionScenario
 * 1 decision variable
 * 1 recourse variable
Solver is default solver

By the law of large numbers, we approach the generalized formulation with increasing sample size. Solving yields:

optimize!(sp, solver = IpoptSolver(print_level=0))

println("Optimal decision: $(optimal_decision(sp))")
println("Optimal value: $(optimal_value(sp))")

Optimal decision: [2.07583]
Optimal value: 4.00553678799426

Now, due to the special implementation of the expected function, it actually holds that the expected value solution solves the generalized problem. Consider:

println("EVP decision: $(EVP_decision(sp, solver = IpoptSolver(print_level=0)))")
println("VSS: $(VSS(sp, solver = IpoptSolver(print_level=0)))")

EVP decision: [2.0]
VSS: 0.005750340653017716

Accordingly, the VSS is small.