Stochastic models

The @stochastic_model command is now introduced in more detail. The discussion will as before revolve around the simple example introduced in the Quick start:

@stochastic_model simple_model begin
    @stage 1 begin
        @decision(simple_model, x₁ >= 40)
        @decision(simple_model, x₂ >= 20)
        @objective(simple_model, Min, 100*x₁ + 150*x₂)
        @constraint(simple_model, x₁ + x₂ <= 120)
    end
    @stage 2 begin
        @known(simple_model, x₁, x₂)
        @uncertain q₁ q₂ d₁ d₂
        @recourse(simple_model, 0 <= y₁ <= d₁)
        @recourse(simple_model, 0 <= y₂ <= d₂)
        @objective(simple_model, Max, q₁*y₁ + q₂*y₂)
        @constraint(simple_model, 6*y₁ + 10*y₂ <= 60*x₁)
        @constraint(simple_model, 8*y₁ + 5*y₂ <= 80*x₂)
    end
end

Note, that the resulting model object is stored in simple_model, and that the same name is used to reference the stochastic program in the @stage blocks. The following anonymous syntax is also supported:

simple_model = @stochastic_model begin
    @stage 1 begin
        @decision(model, x₁ >= 40)
        @decision(model, x₂ >= 20)
        @objective(model, Min, 100*x₁ + 150*x₂)
        @constraint(model, x₁ + x₂ <= 120)
    end
    @stage 2 begin
        @known(model, x₁, x₂)
        @uncertain q₁ q₂ d₁ d₂
        @recourse(model, 0 <= y₁ <= d₁)
        @recourse(model, 0 <= y₂ <= d₂)
        @objective(model, Max, q₁*y₁ + q₂*y₂)
        @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
        @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
    end
end

where the reserved keyword model is used in the @stage blocks.

`@stage` blocks

The body of a @stochastic_model definition consists of a number of @stage blocks, following the syntax:

@stage N begin
    ...
end

Here, N is the stage number and the body is made up of JuMP syntax as well as @parameters, @decision, and @uncertain blocks. At least two stages must be defined and the stages must be defined in consecutive order starting with the first stage. The number of stage blocks included in the @stochastic_model definition determines the number of stages that a stochastic program instantiated from the resulting stochastic model will have.

Note

It is possible to define and instantiate stochastic models with more than two stages. However, most internal tools and solvers only support two-stage models at this point.

`@parameters` blocks

The @parameters blocks are used to introduce deterministic parameters to a @stage block. See for example Stage data. The following:

@parameters a b

makes the constants a and b available as model parameters. This incurs a promise that those parameters will be injected when instantiating the model, and if no default values are available they must be supplied by the user. In other words, if sm is a stochastic model that includes the above @parameters annotation in one of its @stage blocks, then those parameters must be supplied as keyword arguments when instantiating stochastic programs using this model:

instantiate(sm, scenarios, a = 1, b = 2)

Alternatively, default values can be specified directly in the @parameters block:

@parameters begin
    a = 1
    b = 2
end

Values supplied to instantiate are always used, and otherwise the default values are used. The responsibility is on the user to ensure that the supplied parameters support the operations used in the @stage blocks. Parameters can be reused in multiple blocks, but each occurance must be annotated by @parameters in each of the stages.

`@decision` blocks

The @decision blocks are used to annotate linking variables between stages. Their usage is identical syntax-wise to JuMP's @variable macros. Internally, they create specialized JuMP variables with context-dependent behaviour.

`@recourse` blocks

The @recourse macro functions like the @decision macro, but is only used to annotate recourse decisions in the final stage.

`@known` blocks

A @known annotation is used in subsequent stages to bring a decision defined in a previous stage into scope. Any decision defined by @decision inside a @stochastic_model automatically annotates subsequent stages with appropriate @known lines.

The @known block in the simple example above is given by

@known(simple_moel, x₁, x₂)

This states that the second stage of the stochastic model depends on the decisions x₁ and x₂ taken in the previous stage. Note again that this lines is implicitly added by @stochastic_model and is not required.

`@uncertain` blocks

The @uncertain blocks are used to annotate stochastic data in the stochastic model. For flexibility, there are several different ways of doing this. However, an @uncertain annotation is always connected to some AbstractScenario type, as introduced in Scenario data. Note, that a @stage block can only include one @uncertain block. All stochastic information in a given stage must therefore be captured in the @uncertain block of that stage.

The most simple approach is to use leverage the Scenario type. Consider the @uncertain annotation given above:

@uncertain q₁ q₂ d₁ d₂

This will ensure that Scenarios that are expected to have the fields q₁, q₂, d₁ and d₂ are injected when constructing second-stage models. Each such scenario must be supplied or sampled using a supplied sampler object. It is the responsibility of the user to ensure that each supplied or sampled Scenario has the correct fields. For example, the following yields a Scenario compatible with the above @uncertain line:

Scenario(q₁ = 24.0,
         q₂ = 28.0,
         d₁ = 500.0,
         d₂ = 100.0,
         probability = 0.4)

Alternatively, the same scenario is conveniently created using the @scenario macro, matching the syntax of the @uncertain declaration:

@scenario q₁ = 24.0 q₂ = 28.0 d₁ = 500.0 d₂ = 100.0 probability = 0.4

We can also use JuMP's container syntax:

@uncertain ξ[1:5]
@uncertain ξ[i in 1:5]
@uncertain ξ[i in 1:5, i != 3]
@uncertain ξ[i in 1:5, j in 1:5]
@uncertain ξ[i in 1:5, k in [:a,:b,:c]]

and then use the a corresponding formulation in the @scenario macro to generate scenarios:

ξ = @scenario ξ[1:5] = rand(5) probability = rand()
ξ = @scenario ξ[i in 1:5] i * rand() probability = rand()
ξ = @scenario ξ[i in 1:5, i != 3] i * rand() probability = rand()
ξ = @scenario ξ[i in 1:5, j in 1:5] = rand(5,5) probability = rand()
ξ = @scenario ξ[i in 1:5, k in [:a,:b,:c]] = rand(5,5) probability = rand()

Note, that we above sometimes assign the full random vector directly, and sometimes provide an indexed based formula.

As shown in Stochastic data, it is also possible to introduce other scenario types, either using @define_scenario or manally as explained in Custom scenarios and demonstrated in the Continuous scenario distribution example. If we instead define the necessary scenario structure as follows:

@define_scenario SimpleScenario = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end

One can then use:

@uncertain ξ::SimpleScenario

and extract the required fields from ξ which will be of type SimpleScenario after data injection. Again, it is the responsibility of the user to supply scenarios of this type when instantiating the model. For example, the following constructs a SimpleScenario compatible with the above @uncertain line:

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

It is also possible to directly unpack the necessary fields using the following syntactic sugar:

@uncertain q₁ q₂ d₁ d₂ from SimpleScenario

The actual scenario instance can still be annotated and used if necessary:

@uncertain q₁ q₂ d₁ d₂ from ξ::SimpleScenario

Finally, if the @uncertain block is used within a @stochastic_model environment, it is possible to simultaneosly define the underlying scenario type. In other words,

@uncertain ξ::SimpleScenario = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end

@uncertain q₁ q₂ d₁ d₂ from SimpleScenario = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end

and

@uncertain q₁ q₂ d₁ d₂ from ξ::SimpleScenario = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end

are all possible methods of defining and using the SimpleScenario type in a @stage block.

Model instantiation

A model object sm defined using @stochastic_model can be used to instantiate stochastic programs over both finite/infinite sample spaces and discrete/continuous random variables.

If the scenarios are associated with a discrete random variable over a finite sample space, then the corresponding stochastic program is finite and can be instantiated by providing the full list of scenarios:

sp = instantiate(sm, scenarios)

Here, scenarios is a vector of scenarios consistent with the @uncertain annotation used in the second stage of sm. It is the responsibility of the user to ensure that the individual probabilities of the scenarios sum up to one, so that the model is consistent.

If the scenarios are instead associated with a continuous random variable, with finite second moments, over an infinite sample space, then the corresponding stochastic program is not finite and must be approximated. The only supported way of doing so in StochasticPrograms is by using sampled average approximations. A finite stochastic program that approximates the stochastic model is obtained through

sp = instantiate(sm, sampler, n)

where sampler is an AbstractSampler, as outlined in Sampling, and n is the number of samples to include.

Instant models

It is possible to create one-off stochastic programs without needing to first define a model object. To do so, any required scenario data structure must be defined first. Consider:

using StochasticPrograms

@define_scenario SimpleScenario = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end
ξ₁ = SimpleScenario(-24.0, -28.0, 500.0, 100.0, probability = 0.4)
ξ₂ = SimpleScenario(-28.0, -32.0, 300.0, 300.0, probability = 0.6)

SimpleScenario with probability 0.6
  q₁: -28.0
  q₂: -32.0
  d₁: 300.0
  d₂: 300.0

Next, an unmodeled stochastic program can be instantiated using the two created scenarios:

sp = StochasticProgram([ξ₁, ξ₂], Deterministic())

Deferred stochastic program with 2 scenarios

Note that we must provide the instantiation type explicitly as well. A slightly diferrent modeling syntax is now used to define the stage models of sp:

@first_stage sp = begin
    @decision(sp, x₁ >= 40)
    @decision(sp, x₂ >= 20)
    @objective(sp, Min, 100*x₁ + 150*x₂)
    @constraint(sp, x₁ + x₂ <= 120)
end
@second_stage sp = begin
    @known(sp, x₁, x₂)
    @uncertain q₁ q₂ d₁ d₂ from SimpleScenario
    @recourse(sp, 0 <= y₁ <= d₁)
    @recourse(sp, 0 <= y₂ <= d₂)
    @objective(sp, Min, q₁*y₁ + q₂*y₂)
    @constraint(sp, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(sp, 8*y₁ + 5*y₂ <= 80*x₂)
end

Here, @first_stage and @second_stage are just syntactic sugar for @stage 1 and @stage 2. Note, that the model name sp must be used internally in the @stage blocks when referencing the model. This is is the definition syntax used internally by StochasticModel objects when instantiating stochastic programs. Note, that we must explicitly add the @known annotations to the second stage with this approach, while models created using @stochastic_model does this automatically. We can verify that this approach yields the same stochastic program by printing and comparing to the Quick start:

print(sp)

Deterministic equivalent problem
Min 100 x₁ + 150 x₂ - 9.600000000000001 y₁₁ - 11.200000000000001 y₂₁ - 16.8 y₁₂ - 19.2 y₂₂
Subject to
 x₁ ∈ Decisions
 x₂ ∈ Decisions
 y₁₁ ∈ RecourseDecisions
 y₂₁ ∈ RecourseDecisions
 y₁₂ ∈ RecourseDecisions
 y₂₂ ∈ RecourseDecisions
 x₁ ≥ 40.0
 x₂ ≥ 20.0
 y₁₁ ≥ 0.0
 y₂₁ ≥ 0.0
 y₁₂ ≥ 0.0
 y₂₂ ≥ 0.0
 x₁ + x₂ ≤ 120.0
 -60 x₁ + 6 y₁₁ + 10 y₂₁ ≤ 0.0
 -80 x₂ + 8 y₁₁ + 5 y₂₁ ≤ 0.0
 -60 x₁ + 6 y₁₂ + 10 y₂₂ ≤ 0.0
 -80 x₂ + 8 y₁₂ + 5 y₂₂ ≤ 0.0
 y₁₁ ≤ 500.0
 y₂₁ ≤ 100.0
 y₁₂ ≤ 300.0
 y₂₂ ≤ 300.0
Solver name: No optimizer attached.

As a side note, it is possible to run stage definition macros on programs with existing models. This overwrites the previous model and upon regeneration all internal problems. For example, the following increases the lower bound on the second stage variables to 2:

@stage 2 sp = begin
    @known(sp, x₁, x₂)
    @uncertain q₁ q₂ d₁ d₂ from SimpleScenario
    @recourse(sp, 2 <= y₁ <= d₁)
    @recourse(sp, 2 <= y₂ <= d₂)
    @objective(sp, Min, q₁*y₁ + q₂*y₂)
    @constraint(sp, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(sp, 8*y₁ + 5*y₂ <= 80*x₂)
end

generate!(sp)

print(sp)

Deterministic equivalent problem
Min 100 x₁ + 150 x₂ - 9.600000000000001 y₁₁ - 11.200000000000001 y₂₁ - 16.8 y₁₂ - 19.2 y₂₂
Subject to
 x₁ ∈ Decisions
 x₂ ∈ Decisions
 y₁₁ ∈ RecourseDecisions
 y₂₁ ∈ RecourseDecisions
 y₁₂ ∈ RecourseDecisions
 y₂₂ ∈ RecourseDecisions
 x₁ ≥ 40.0
 x₂ ≥ 20.0
 y₁₁ ≥ 2.0
 y₂₁ ≥ 2.0
 y₁₂ ≥ 2.0
 y₂₂ ≥ 2.0
 x₁ + x₂ ≤ 120.0
 -60 x₁ + 6 y₁₁ + 10 y₂₁ ≤ 0.0
 -80 x₂ + 8 y₁₁ + 5 y₂₁ ≤ 0.0
 -60 x₁ + 6 y₁₂ + 10 y₂₂ ≤ 0.0
 -80 x₂ + 8 y₁₂ + 5 y₂₂ ≤ 0.0
 y₁₁ ≤ 500.0
 y₂₁ ≤ 100.0
 y₁₂ ≤ 300.0
 y₂₂ ≤ 300.0
Solver name: No optimizer attached.

It is of course also possible to do this on programs instantiated from a StochasticModel.

SMPS

StochasticPrograms also support reading models specified in the SMPS format. Specifically, SMPS definitions with uncertain data of types INDEP or BLOCKS are supported. We show how the simple example can be specified in SMPS. An SMPS definition consist of the following files:

problem.smps
problem.tim
problem.cor
problem.sto

Here, problem.smps is just an empty file that shares the name with the others to simplify IO commands. The problem.tim file specifies the stage structure of the stochastic program. An example for the simple problem is given below.

TIME          SIMPLE
PERIODS
    X1        BOUND                    STAGE1
    Y1        LINK1                    STAGE2
ENDATA

Row and column delimeters are given for each stage. The problem.cor file specifies the optimization structure of the problem in MPS format. An example for the simple problem as follows.

NAME          SIMPLE
ROWS
 N  OBJ
 L  BOUND
 L  LINK1
 L  LINK2
 L  Y1UP
 L  Y2UP
COLUMNS
    X1        OBJ       100.0          BOUND     1.0
    X2        OBJ       150.0          BOUND     1.0
    X1        LINK1     -60.0
    X2        LINK2     -80.0
    Y1        OBJ       26.0           Y1UP      1.0
    Y2        OBJ       30.0           Y2UP      1.0
    Y1        LINK1     6.0            LINK2     8.0
    Y2        LINK1     10.0           LINK2     5.0
RHS
    RHS       BOUND     120.0
    RHS       LINK1     0.0
    RHS       LINK2     0.0
    RHS       Y1UP      400.0
    RHS       Y2UP      200.0
BOUNDS
 LO X1LIM     X1        40.0
 LO X2LIM     X2        20.0
 LO Y1LIM     Y1        0.0
 LO Y2LIM     Y2        0.0
ENDATA

Finally, the problem.sto file specifies the uncertain data. We use the BLOCKS format to specify the simple scenarios.

STOCH         SIMPLE
BLOCKS        DISCRETE
 BL BLOCK1    STAGE2    0.4
    Y1        OBJ       -24.0
    Y2        OBJ       -28.0
    RHS       Y1UP      500.0
    RHS       Y2UP      100.0
 BL BLOCK1    STAGE2    0.6
    Y1        OBJ       -28.0
    Y2        OBJ       -32.0
    RHS       Y1UP      300.0
    RHS       Y2UP      300.0
ENDATA

We specify the two scenarios, giving the value in each corresponding coordinate in the problem.cor file. Now, we can read this model into Julia in several ways, assuming all files are in the same folder. First, consider

model = read("problem.smps", StochasticModel)

Two-Stage Stochastic Model

minimize f₀(x) + 𝔼[f(x,ξ)]
  x∈𝒳

where

f(x,ξ) = min  f(y; x, ξ)
              y ∈ 𝒴 (x, ξ)

This returns a StochasticModel object that can be used as usual, assuming it is instantiated with the scenarios of the special SMPSScenario type. To that end, we can read a specialized sampler object for the specified SMPS model:

sampler = read("problem.smps", SMPSSampler)
sampler()

SMPSScenario with probability 1.0 and underlying data:

Δq = [0.0, 0.0, -54.0, -62.0]
ΔT = 0×2 SparseArrays.SparseMatrixCSC{Float64,Int64} with 0 stored entries
ΔW = 0×2 SparseArrays.SparseMatrixCSC{Float64,Int64} with 0 stored entries
Δh = Float64[]
ΔC = 4×4 SparseArrays.SparseMatrixCSC{Float64,Int64} with 0 stored entries
Δd₁ = [0.0, 0.0, 0.0, 0.0]
Δd₂ = [0.0, 0.0, -100.0, 100.0]

We can now instantiate a specific instance of the read model:

sp = instantiate(model, sampler, 2)

Stochastic program with:
 * 2 decision variables
 * 2 recourse variables
 * 2 scenarios of type SMPSScenario
Structure: Deterministic equivalent
Solver name: No optimizer attached.

The same result, up to sampling, can be achieved directly through

sp = read("problem.smps", StochasticProgram, num_scenarios = 2)

Stochastic program with:
 * 2 decision variables
 * 2 recourse variables
 * 2 scenarios of type SMPSScenario
Structure: Deterministic equivalent
Solver name: No optimizer attached.

This read variant takes the same keyword arguments as instantiate. Because the specified scenario structure in BLOCKS or INDEP format has finite support, it is possible to read the stochastic program corresponding to the full support by not specifying num_scenarios:

sp = read("problem.smps", StochasticProgram)

Stochastic program with:
 * 2 decision variables
 * 2 recourse variables
 * 2 scenarios of type SMPSScenario
Structure: Deterministic equivalent
Solver name: No optimizer attached.

In this case, the full support correspond exactly to the simple model we have considered before:

print(sp)

Deterministic equivalent problem
Min 100 x[1] + 150 x[2] - 16.8 y₂[1] - 19.2 y₂[2] - 9.600000000000001 y₁[1] - 11.200000000000001 y₁[2]
Subject to
 y₁[1] ≥ 0.0
 y₁[2] ≥ 0.0
 y₂[1] ≥ 0.0
 y₂[2] ≥ 0.0
 [x[1], x[2]] ∈ Decisions
 x[1] ≥ 40.0
 x[2] ≥ 20.0
 [x[1] + x[2] - 120] ∈ MathOptInterface.Nonpositives(1)
 [-60 x[1] + 6 y₁[1] + 10 y₁[2], -80 x[2] + 8 y₁[1] + 5 y₁[2], y₁[1] - 500, y₁[2] - 100] ∈ MathOptInterface.Nonpositives(4)
 [-60 x[1] + 6 y₂[1] + 10 y₂[2], -80 x[2] + 8 y₂[1] + 5 y₂[2], y₂[1] - 300, y₂[2] - 300] ∈ MathOptInterface.Nonpositives(4)
Solver name: No optimizer attached.

A warning is issued if the full support contains more than 1e5 scenarios.

Stochastic models

@stage blocks

@parameters blocks

@decision blocks

@recourse blocks

@known blocks

@uncertain blocks

Model instantiation

Instant models

SMPS

`@stage` blocks

`@parameters` blocks

`@decision` blocks

`@recourse` blocks

`@known` blocks

`@uncertain` blocks