Suppose I have a two stage problem where I expect the user to provide (1) a first-stage model written in JuMP and (2) a function to evaluate the recourse. If the user’s recourse function evaluates subproblems also written in JuMP, then I would like to formulate the extensive form problem by copying the first-stage model and then copying variables/constraints/objectives from the second stage into the first stage.
Right now, I am able to copy objectives of models 1...n into a single model. But I am having trouble doing the same for the constraints. My attempt (below; initially here) was to collect the expressions of variables/objectives/constraints and try to remap them.
I’ve been able to extract the objectives as expressions and remap variables from the collection of models to a new monolithic model. However, I have trouble on constraints. E.g., if a model i has @variable(modeli, y[1] >=0) and @constraint(modeli, y[1] >= 1), my method is not sure how to handle this.
Is there a nice way of doing this for constraints? Or should I consider a different API for solving the extensive form model?
MWE:
using JuMP, OrderedCollections
#=
start utilities --------------------------------------------------------------------------------------------------------------
=#
function substitute_variables(expr, mapping::Dict{JuMP.VariableRef, JuMP.VariableRef})
if expr isa JuMP.QuadExpr
# Handle Quadratic Expressions
new_quad_terms = OrderedCollections.OrderedDict(
UnorderedPair(get(mapping, pair.a, pair.a), get(mapping, pair.b, pair.b)) => coef
for (pair, coef) in expr.terms
)
# Recursively substitute the affine part
new_aff_expr = substitute_variables(expr.aff, mapping)
return JuMP.QuadExpr(new_aff_expr, new_quad_terms)
elseif expr isa JuMP.AffExpr
# Handle Affine Expressions
new_terms = OrderedCollections.OrderedDict(
get(mapping, var, var) => coef for (var, coef) in expr.terms
)
return JuMP.AffExpr(expr.constant, new_terms...)
elseif expr isa JuMP.VariableRef
# Handle Single Variables
return get(mapping, expr, expr)
else
# nonlinear expressions...see: https://jump.dev/JuMP.jl/stable/manual/nonlinear/#Nonlinear-expressions-in-detail
error("Unsupported expression type: $(typeof(expr))")
end
end
#=
end utilities --------------------------------------------------------------------------------------------------------------
=#
# Step 1: Define individual models
n = 3 # Number of individual models
models = []
for i in 1:n
model = Model()
@variable(model, y[1:2] >= 0) # Each model has a variable array `y` of size 2
@objective(model, Min, y[1]^2 + i * y[2]) # Example objective
@constraint(model, y[1] + i <= 10) # Example constraint
@constraint(model, y[2] + i >= 10) # Example constraint
push!(models, model)
end
# Step 2: Create the big model
big_model = Model()
# Step 3: Create variables for the big model and build a variable mapping
variable_mapping = Dict{JuMP.VariableRef, JuMP.VariableRef}()
for (model_idx, model) in enumerate(models)
# Iterate over all variables in the model
for var in all_variables(model)
# Get the variable's symbolic name
var_name = Symbol(var) # e.g., "y[1]", "y[2]"
# Create a unique name for the big model variable
unique_var_name = Symbol("m$(model_idx)_$(var_name)")
# Get the variable's properties
lb = JuMP.has_lower_bound(var) ? JuMP.lower_bound(var) : -Inf
ub = JuMP.has_upper_bound(var) ? JuMP.upper_bound(var) : Inf
is_bin = JuMP.is_binary(var)
is_int = JuMP.is_integer(var)
# Create the variable in the big model with the same properties
if JuMP.has_lower_bound(var) && JuMP.has_upper_bound(var)
new_var = @variable(big_model, lower_bound=lb, upper_bound=ub, base_name=string(unique_var_name))
elseif JuMP.has_lower_bound(var) && !JuMP.has_upper_bound(var)
new_var = @variable(big_model, lower_bound=lb, base_name=string(unique_var_name))
elseif !JuMP.has_lower_bound(var) && JuMP.has_upper_bound(var)
new_var = @variable(big_model, upper_bound=ub, base_name=string(unique_var_name))
else
new_var = @variable(big_model, base_name=string(unique_var_name))
end
# Map the small model variable to the newly created big model variable
variable_mapping[var] = new_var
end
end
# Verify the variable mapping
println("Variable Mapping:")
for (small_var, big_var) in variable_mapping
println("Small model variable $small_var -> Big model variable $big_var")
end
# -------------------------------------------------------------------------------------------------------------------------------
# Step 4: Combine objectives and constraints from all models
combined_objective = 0.0 # Initialize the combined objective
for (model_idx, model) in enumerate(models)
# Extract the objective function
obj_expr = JuMP.objective_function(model)
# Perform substitution for variables using the variable_mapping
substituted_obj_expr = substitute_variables(obj_expr, variable_mapping)
# Add the substituted objective to the combined objective
combined_objective += substituted_obj_expr
end
# Add the combined objective to the big model
@objective(big_model, Min, combined_objective)
# Step 5: Extract and add constraints from individual models to the big model
for (model_idx, model) in enumerate(models)
println("model=$model_idx")
for (con_idx, con) in enumerate(JuMP.all_constraints(model, include_variable_in_set_constraints=true))
println("con=$con, $con_idx")
# Get the full constraint object
constraint_obj = JuMP.constraint_object(con)
# Access the left-hand side (expression) and right-hand side (set)
con_expr = constraint_obj.func # The function or expression on the LHS
@show con_expr
con_set = constraint_obj.set # The set defining the constraint (e.g., `<=`, `>=`, `==`)
@show con_set
# Substitute variables in the constraint expression
substituted_con_expr = substitute_variables(con_expr, variable_mapping)
# Add the substituted constraint to the big model
@constraint(big_model, substituted_con_expr in con_set)
end
end
# Verify the big model
println("Big Model Objective:")
println(combined_objective)
println("Big Model Constraints:")
for con in JuMP.all_constraints(big_model)
println(con)
end
gives an error
model=1
con=y[2] ≥ 9, 1
con_expr = y[2]
con_set = MathOptInterface.GreaterThan{Float64}(9.0)
con=y[1] ≤ 9, 2
con_expr = y[1]
con_set = MathOptInterface.LessThan{Float64}(9.0)
con=y[1] ≥ 0, 3
con_expr = y[1]
con_set = MathOptInterface.GreaterThan{Float64}(0.0)
ERROR: MathOptInterface.LowerBoundAlreadySet{MathOptInterface.GreaterThan{Float64}, MathOptInterface.GreaterThan{Float64}}: Cannot add `VariableIndex`-in-`MathOptInterface.GreaterThan{Float64}` constraint for variable MOI.VariableIndex(1) as a `VariableIndex`-in-`MathOptInterface.GreaterThan{Float64}` constraint was already set for this variable and both constraints set a lower bound.
Stacktrace:
[1] _throw_if_lower_bound_set_inner(variable::MathOptInterface.VariableIndex, S2::Type, mask::UInt16, T::Type)
@ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/variables_container.jl:120
[2] _throw_if_lower_bound_set
@ ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/variables_container.jl:131 [inlined]
[3] add_constraint
@ ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/variables_container.jl:267 [inlined]
[4] add_constraint
@ ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/model.jl:348 [inlined]
[5] add_constraint
@ ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/universalfallback.jl:835 [inlined]
[6] add_constraint(m::MathOptInterface.Utilities.CachingOptimizer{…}, func::MathOptInterface.VariableIndex, set::MathOptInterface.GreaterThan{…})
@ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/cachingoptimizer.jl:551
[7] _moi_add_constraint(model::MathOptInterface.Utilities.CachingOptimizer{…}, f::MathOptInterface.VariableIndex, s::MathOptInterface.GreaterThan{…})
@ JuMP ~/.julia/packages/JuMP/V9nZm/src/constraints.jl:1013
[8] add_constraint(model::Model, con::ScalarConstraint{VariableRef, MathOptInterface.GreaterThan{Float64}}, name::String)
@ JuMP ~/.julia/packages/JuMP/V9nZm/src/constraints.jl:1036
[9] macro expansion
@ ~/.julia/packages/JuMP/V9nZm/src/macros/@constraint.jl:173 [inlined]
[10] macro expansion
@ ~/.julia/packages/JuMP/V9nZm/src/macros.jl:393 [inlined]
[11] top-level scope
@ ./REPL[165]:18
Some type information was truncated. Use `show(err)` to see complete types.
