Decision API
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(sp, x >= 2)
@variable(sp, w)
@objective(sp, Min, x)
end
@second_stage sp = begin
@known(sp, x)
@uncertain a
@recourse(sp, y >= 1/a)
@objective(sp, Max, y)
@constraint(sp, 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₁ : -x + y₁ ≤ 0.0
con₂ : -x + 2 y₂ ≤ 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 Vector{DecisionVariable}:
x
and
all_decision_variables(sp, 2)
1-element Vector{DecisionVariable}:
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) = 1.0
y₂
y has lower bound in scenario 2: true
y has upper bound in scenario 2: false
lower_bound(y, 2) = 0.5
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 StageDecomposition
or ScenarioDecomposition
structure. Auxiliary 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{Type, Type}[(DecisionRef, MathOptInterface.GreaterThan{Float64})]
println(list_of_constraint_types(sp, 2))
Tuple{Type, Type}[(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 Vector{SPConstraintRef{MathOptInterface.ConstraintIndex{AffineDecisionFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:
con : -x + y ≤ 0.0
The scenario-dependent constraint in stage 2 can also be accessed through
con = sp[2,:con];
println(con)
con : -x + y ≤ 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 instead instantiated with StageDecomposition
or ScenarioDecomposition
.
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))")
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₁ : -x + y₁ ≤ 0.0
con₂ : -x + 2 y₂ ≤ 0.0
Solver name: GLPK
This is more apparent in a stage-decomposition structure:
vertical_sp = StochasticProgram([ξ₁,ξ₂], StageDecomposition())
@first_stage vertical_sp = begin
@decision(vertical_sp, x >= 2)
@variable(vertical_sp, w)
@objective(vertical_sp, Min, x)
end
@second_stage vertical_sp = begin
@known(vertical_sp, x)
@uncertain a
@recourse(vertical_sp, y >= 1/a)
@objective(vertical_sp, Max, y)
@constraint(vertical_sp, 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 : -x + y ≤ 0.0
Subproblem 2 (p = 0.50):
Max y
Subject to
x ∈ Known(value = 2.0)
y ∈ RecourseDecisions
y ≥ 0.5
con : -x + 2 y ≤ 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.]))")
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: 3.25
Optimal decision: [3.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.5
Objective value in scenario 2: 1.5
Optimal recourse in scenario 2: [1.5]
3.25
Equivalent decision evaluation: 3.25