Multi-stage decisions

StochasticPrograms provides an extended version of JuMPs API for decision variables/constraints/objectives with stage and scenario dependence. Consider the following problem:

using StochasticPrograms
using GLPK

ξ₁ = @scenario a = 1 probability = 0.5
ξ₂ = @scenario a = 2 probability = 0.5

sp = StochasticProgram([ξ₁,ξ₂], Deterministic())
@first_stage sp = begin
    @decision(model, x >= 2)
    @variable(model, w)
    @objective(model, Min, x)
end
@second_stage sp = begin
    @known x
    @uncertain a
    @recourse(model, y >= 1/a)
    @objective(model, Max, y)
    @constraint(model, con, a*y <= x)
end

A first-stage decision and a second-stage recourse decision has been identified. We include a printout of the stochastic program for reference:

print(sp)

Deterministic equivalent problem
Min x - 0.5 y₁ - 0.5 y₂
Subject to
 x ∈ Decisions
 y₁ ∈ RecourseDecisions
 y₂ ∈ RecourseDecisions
 x ≥ 2.0
 y₁ ≥ 1.0
 y₂ ≥ 0.5
 con : y₁ - x ≤ 0.0
 con : 2 y₂ - x ≤ 0.0
Solver name: No optimizer attached.

Decision variables

All variables annotated with either @decision or @recourse become available through the API. We can query all such variables:

all_decision_variables(sp)

(DecisionVariable[x], DecisionVariable[y])

stage-wise lists can also be obtained:

all_decision_variables(sp, 1)

1-element Array{DecisionVariable,1}:
 x

and

all_decision_variables(sp, 2)

1-element Array{DecisionVariable,1}:
 y

JuMPs [] syntax is available as well, with the addition that the stage must be provided as well:

x = sp[1,:x]
println(x)
println(typeof(x))

x
DecisionVariable

The return type is DecisionVariable a specialized AbstractVariableRef. For first-stage variables, the syntax is unchanged from JuMP:

println(name(x))
println("x has lower bound: $(has_lower_bound(x))")
println("x has upper bound: $(has_upper_bound(x))")
println("lower_bound(x) = $(lower_bound(x))")

x
x has lower bound: true
x has upper bound: false
lower_bound(x) = 2.0

If we instead query for the recourse decision $y$:

y = sp[2,:y]
println(y)
println(typeof(y))

y
DecisionVariable

The same getters will error:

lower_bound(y)

ERROR: y is scenario dependent, consider `lower_bound(dvar, scenario_index)`.

As indicated by the error, $y$ is scenario-dependent so a scenario_index must be provided as well:

println(name(y, 1))
println("y has lower bound in scenario 1: $(has_lower_bound(y, 1))")
println("y has upper bound in scenario 1: $(has_upper_bound(y, 1))")
println("lower_bound(y, 1) = lower_bound(y, 1)")

println(name(y, 2))
println("y has lower bound in scenario 2: $(has_lower_bound(y, 2))")
println("y has upper bound in scenario 2: $(has_upper_bound(y, 2))")
println("lower_bound(y, 2) = lower_bound(y, 2)")

y₁
y has lower bound in scenario 1: true
y has upper bound in scenario 1: false
lower_bound(y, 1) = lower_bound(y, 1)
y₂
y has lower bound in scenario 2: true
y has upper bound in scenario 2: false
lower_bound(y, 2) = lower_bound(y, 2)

The lower bound of $y$ is as expected different in the two scenarios. Some attributes, such as the variable name, are structure dependent and may vary in a Vertical or Horizontal structure. Auxilliary variables created with the standard @variable are not available through this API. To access them, either annotate them with @decision (or @recourse in the final stage), or access the relevant JuMP subproblem and query the variable as usual. For example:

w = DEP(sp)[:w]
println(w)
println(typeof(w))

w
VariableRef

Decision constraints

Constraints that include variables annotated with either @decision or @recourse can also be accessed in the extended API. Stage-wise list of all such constraints can be obtained:

println(list_of_constraint_types(sp, 1))

Tuple{DataType,DataType}[(DecisionRef, MathOptInterface.GreaterThan{Float64})]

println(list_of_constraint_types(sp, 2))

Tuple{DataType,DataType}[(DecisionRef, MathOptInterface.GreaterThan{Float64}), (DecisionAffExpr{Float64}, MathOptInterface.LessThan{Float64})]

and type-sorted constraints can be obtained through a stage-dependent variant of all_constraints:

all_constraints(sp, 2, DecisionAffExpr{Float64}, MOI.LessThan{Float64});

1-element Array{SPConstraintRef{MathOptInterface.ConstraintIndex{AffineDecisionFunction{Float64},MathOptInterface.LessThan{Float64}},ScalarShape},1}:
 con : w - x ≤ 0.0

The scenario-dependent constraint in stage 2 can also be accessed through

con = sp[2,:con];
println(con)

con : w - x ≤ 0.0

This returns an SPConstraintRef, similar in function to DecisionVariable. The constraint originates from stage-two, so most attributes are scenario-dependent:

println(name(con, 1))
println("RHS of con in scenario 1 = $(normalized_rhs(con, 1))")
println("Coefficient of x in scenario 1 = $(normalized_coefficient(con, x, 1))")
println("Coefficient of y in scenario 1 = $(normalized_coefficient(con, y, 1))")

println(name(con, 2))
println("RHS of con in scenario 2 = $(normalized_rhs(con, 2))")
println("Coefficient of x in scenario 2 = $(normalized_coefficient(con, x, 2))")
println("Coefficient of y in scenario 2 = $(normalized_coefficient(con, y, 2))")

con
RHS of con in scenario 1 = 0.0
Coefficient of x in scenario 1 = -1.0
Coefficient of y in scenario 1 = 1.0
con
RHS of con in scenario 2 = 0.0
Coefficient of x in scenario 2 = -1.0
Coefficient of y in scenario 2 = 2.0

Decision objectives

The objective function of a stochastic program can be obtained in full or in stage and scenario-dependent chunks:

println("Objective in stage 1: $(objective_function(sp, 1))")
println("Objective in stage 2, scenario 1: $(objective_function(sp, 2, 1))")
println("Objective in stage 2, scenario 2: $(objective_function(sp, 2, 2))")
println("Full objective: $(objective_function(sp))")

Objective in stage 1: x
Objective in stage 2, scenario 1: y₁
Objective in stage 2, scenario 2: y₂
Full objective: x - 0.5 y₁ - 0.5 y₂

and can be modified accordingly:

set_objective_coefficient(sp, y, 2, 1, 2.);
println("Objective in stage 2, scenario 1: $(objective_function(sp, 2, 1))")
println("Full objective: $(objective_function(sp))")
set_objective_coefficient(sp, y, 2, 1, 1.);

Objective in stage 2, scenario 1: 2 y₁
Full objective: x + y₁ - 0.5 y₂

The stochastic program objective is structure dependent and will appear different if the stochastic program is instantiated with Vertical or Horizontal instead.

Solved problem

After optimizing the stochastic program, attributes for which is_set_by_optimize is true can be accessed using the same scenario-dependent syntax:

set_optimizer(sp, GLPK.Optimizer)

optimize!(sp)

# Main result
println("Termination status: $(termination_status(sp))")
println("Objective value: $(objective_value(sp))")
println("Optimal decision: $(optimal_decision(sp))")

# First stage
println("value(x) = $(value(x))")
println("reduced_cost(x) = $(reduced_cost(x))")

# Scenario 1
# Second stage
println("value(y, 1) = $(value(y, 1))")
println("reduced_cost(y, 1) = $(reduced_cost(y, 1))")
println("dual(con, 1) = $(dual(con, 1))")
println("Objective value in scenario 1: $(objective_value(sp, 1))")
println("Optimal recourse in scenario 1: $(optimal_recourse_decision(sp, 1))")

# Scenario 2
println("value(y, 2) = $(value(y, 2))")
println("reduced_cost(y, 2) = $(reduced_cost(y, 2))")
println("dual(con, 2) = $(dual(con, 2))")
println("Objective value in scenario 2: $(objective_value(sp, 2))")
println("Optimal recourse in scenario 2: $(optimal_recourse_decision(sp, 2))")

Termination status: OPTIMAL
Objective value: 2.0
Optimal decision: [2.0]
value(x) = 2.0
reduced_cost(x) = 0.75
value(y, 1) = 1.0
reduced_cost(y, 1) = 0.5
dual(con, 1) = 0.0
Objective value in scenario 1: 2.0
Optimal recourse in scenario 1: [1.0]
value(y, 2) = 1.0
reduced_cost(y, 2) = 0.0
dual(con, 2) = -0.25
Objective value in scenario 2: 2.0
Optimal recourse in scenario 2: [1.0]

As mentioned in the Quick start, decision evaluation can be performed manually through the decision API. Consider:

fix(x, 3.);

This not only fixes $x$ in the first-stage, but also in all occurances in subsequent stages:

print(sp)

Deterministic equivalent problem
Min x + 0.5 y₁ - 0.5 y₂
Subject to
 x ∈ Decisions(value = 3.0)
 y₁ ∈ RecourseDecisions
 y₂ ∈ RecourseDecisions
 x ≥ 2.0
 y₁ ≥ 1.0
 y₂ ≥ 0.5
 con : y₁ - x ≤ 0.0
 con : 2 y₂ - x ≤ 0.0
Solver name: GLPK

This is more apparent in a vertical structure:

vertical_sp = StochasticProgram([ξ₁,ξ₂], Vertical())
@first_stage vertical_sp = begin
    @decision(model, x >= 2)
    @variable(model, w)
    @objective(model, Min, x)
end
@second_stage vertical_sp = begin
    @known x
    @uncertain a
    @recourse(model, y >= 1/a)
    @objective(model, Max, y)
    @constraint(model, con, a*y <= x)
end

x = vertical_sp[1,:x]
fix(x, 3.)

print(vertical_sp)

First-stage
==============
Min x
Subject to
 x ∈ Decisions(value = 3.0)
 x ≥ 2.0

Second-stage
==============
Subproblem 1 (p = 0.50):
Max y
Subject to
 x ∈ Known(value = 2.0)
 y ∈ RecourseDecisions
 y ≥ 1.0
 con : y - x ≤ 0.0

Subproblem 2 (p = 0.50):
Max y
Subject to
 x ∈ Known(value = 2.0)
 y ∈ RecourseDecisions
 y ≥ 0.5
 con : 2 y - x ≤ 0.0

Solver name: No optimizer attached.

We resolve the problem to verify:

optimize!(sp)

# Main result
println("Termination status: $(termination_status(sp))")
println("Objective value: $(objective_value(sp))")
println("Optimal decision: $(optimal_decision(sp))")

# First stage
println("value(x) = $(value(x))")

# Scenario 1
# Second stage
println("value(y, 1) = $(value(y, 1))")
println("reduced_cost(y, 1) = $(reduced_cost(y, 1))")
println("dual(con, 1) = $(dual(con, 1))")
println("Objective value in scenario 1: $(objective_value(sp, 1))")
println("Optimal recourse in scenario 1: $(optimal_recourse_decision(sp, 1))")

# Scenario 2
println("value(y, 2) = $(value(y, 2))")
println("reduced_cost(y, 2) = $(reduced_cost(y, 2))")
println("dual(con, 2) = $(dual(con, 2))")
println("Objective value in scenario 2: $(objective_value(sp, 2))")
println("Optimal recourse in scenario 2: $(optimal_recourse_decision(sp, 2))")

# Evaluating x = 3 should give the same answer:
println("Equivalent decision evaluation: $(evaluate_decision(sp, [3.]))")

Termination status: OPTIMAL
Objective value: 2.75
Optimal decision: [3.0]
value(x) = 3.0
value(y, 1) = 1.0
reduced_cost(y, 1) = 0.5
dual(con, 1) = 0.0
Objective value in scenario 1: 3.0
Optimal recourse in scenario 1: [1.0]
value(y, 2) = 1.5
reduced_cost(y, 2) = 0.0
dual(con, 2) = -0.25
Objective value in scenario 2: 2.75
Optimal recourse in scenario 2: [1.5]
Equivalent decision evaluation: 2.75

We can also fix the value of $y$ in a specific scenario:

fix(y, 1, 2.)
optimize!(sp)

# Main result
println("Termination status: $(termination_status(sp))")
println("Objective value: $(objective_value(sp))")
println("Optimal decision: $(optimal_decision(sp))")

# First stage
println("value(x) = $(value(x))")

# Scenario 1
# Second stage
println("value(y, 1) = $(value(y, 1))")
println("Objective value in scenario 1: $(objective_value(sp, 1))")
println("Optimal recourse in scenario 1: $(optimal_recourse_decision(sp, 1))")

# Scenario 2
println("value(y, 2) = $(value(y, 2))")
println("Objective value in scenario 2: $(objective_value(sp, 2))")
println("Optimal recourse in scenario 2: $(optimal_recourse_decision(sp, 2))")

# Evaluating x = 3 should give the same answer:
println(evaluate_decision(sp, [3.]))

# Evaluating x = 3 should give the same answer:
println("Equivalent decision evaluation: $(evaluate_decision(sp, [3.]))")

Termination status: OPTIMAL
Objective value: 2.5
Optimal decision: [2.0]
value(x) = 3.0
value(y, 1) = 2.0
Objective value in scenario 1: 2.0
Optimal recourse in scenario 1: [2.0]
value(y, 2) = 1.0
Objective value in scenario 2: 2.5
Optimal recourse in scenario 2: [1.0]
3.25
Equivalent decision evaluation: 3.25