Have an strange error creating optimization subproblems on Julia

[benjamin.nunez]

Hello! I’m experiencing a rather strange problem in Julia when creating optimization subproblems.

Basically I have a function called “build_subproblem(s)” which depends on a scenario “s” and formulates an optimization problem using JuMP and specific data related to that scenario “s”.

When I want to create the subproblems for all of the “s” scenarios, I typically write “for s in 1:scenarios” and store the subproblems either in a dictionary or in an array of type Model.

The problem comes once I have already created the subproblems in said cycle, I have another function called “solve_subproblem(Model::)” which receives the subproblems and optimizes them, however I must use it within the same cycle in which I create the subproblems , if I use the function in another line of code or in another block of the notebook I get the error “OptimizeNotCalled()”, which does not make much sense since within the function “solve_subproblem(Model::)” precisely the subproblem is optimized . This means that every time I want to work on a subproblem I have to create it, while the ideal would be to create all the subproblems once and then do whatever I need with them later.

[mkitti]

Could you please provide a minimum working example that allows us to reproduce your problem?

Please read: make it easier to help you

Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …

Perhaps the way you are creating the subproblems is not as independent as you think.

[odow]

Hi @benjamin.nunez, welcome to the forum

As @mkitti says, it would be helpful if you could provide a minimal reproducible example.

My guess is that you’re querying the value of a variable that belongs to a model other than the one that you just solved. But that’s just a guess. If you can provide an example, we should be able to help further.

It’s likely just a small bug. JuMP is well-suited to iterative algorithms where you solve a bunch of scenario-dependent subproblems, see, for example, GitHub - odow/SDDP.jl: Stochastic Dual Dynamic Programming in Julia.

[benjamin.nunez]

Hello @odow ! I appreciate your response and your interest in helping me.

I will share with you the fragments of my code that involve the error I mentioned. In particular, i am using the following data from a simple example of stochastic two-stage linear planning: (just to give you some context) This involves an optimal planning problem for generation expansion, which means determining the optimal investment in generation. There are various types of candidate generators available to meet demand over several periods. The first period has random demand. The investment decision must be consistent across all periods. Additionally, there is a budget constraint on total investment and a minimum power requirement for installation.

# Data:

# Number of scenarios of the problem
scenarios = 3
# Probability of each scenario
probabiliity = [0.2, 0.5, 0.3]

# Demand for each period j under a scenario s
demand = []
push!(demand, [3.0, 3.0, 2.0])
push!(demand, [5.0, 3.0, 2.0])
push!(demand, [7.0, 3.0, 2.0])

# Number of generators 
generators = 4 
# Number of periods
periods = 3 

# Fixed cost of installing generator i
fixed_cost = [10, 7, 16, 6]

# Variable cost of production of generator i in period j
variable_cost =[40 24 4;
                45 27 4.5;
                32 19.2 3.2;
                55 33 5.5]
                
# Minimum power to be installed            
minimum_power = 12

# Budget of the decision maker
budget = 120

Then, the following function takes an input scenario index "s" and is responsible for creating the subproblem (or the deterministic version of the stochastic model for that scenario).

# Dictionary that will store the first-stage variables of the stochastic problem
first_stage_variables = Dict()

# Function that will build a subproblem depending of the data of the scenario "s"
function Build_Subproblem(s::Int)

    # Attributes of the solver (Gurobi)
    subproblem = Model(Gurobi.Optimizer)
    set_optimizer_attribute(subproblem, "OutputFlag", 0)
    set_optimizer_attribute(subproblem, "LogFile", "")
    set_optimizer_attribute(subproblem, "LogToConsole", 0)

    # Decision variables:
    @variable(subproblem, x[1:generators] >= 0)
    @variable(subproblem, y[1:generators, 1:periods] >= 0)

    # We add the first-stage variables to the dictionary:
    first_stage_variables[1] = x
    
    # Subproblem constraints for any given scenario:
    @constraint(subproblem, sum(x[i] for i in 1:generators) >= minimum_power)

    @constraint(subproblem, sum(fixed_cost[i] * x[i] for i in 1:generators) <= budget)
    
    for i in 1:generators
        for j in 1:periods
            @constraint(subproblem, y[i, j] <= x[i])
        end
    end

    for j in 1:periods
        @constraint(subproblem, sum(y[i, j] for i in 1:generators) >= demand[s][j])
    end

    # Objective function of the subproblem for any given scenario:
    @objective(subproblem, Min, sum(fixed_cost[i] * x[i] for i in 1:generators) + sum(sum(variable_cost[i, j] * y[i, j] for i in 1:generators) for j in 1:periods))
    
    # The function returns the subproblem (an optimization model created with JuMP):
    return subproblem
end

And finally, this function takes as input an optimization model (in particular, the subproblems we have created) and is responsible for solving it and delivering the solution for the decision variables of the first stage for the corresponding scenario.

function Solve_Subproblem(subproblem::Model)

    # First we optimize the subproblem:
    optimize!(subproblem)
    
    # An array is created to store the solutions of the first-stage variables:
    first_stage_solutions = []

    # Let's go through all the variables of the first stage:
    for i in 1:length(first_stage_variables)
        # We call the variables saved in the dictionary:
        variable = first_stage_variables[i]
        # We save the solution of the first stage variable:
        push!(first_stage_solutions, JuMP.value.(variable))
    end

    # The solutions of all the variables of the first stage of the subproblem are delivered:
    return first_stage_solutions
end

Now is when you can see the error I mentioned. I would like to be able to create all the subproblems only once and then work with them, for example:

# We create all the subproblems:
subproblems = Model[]
for s in 1:scenarios
    push!(subproblems, Build_Subproblem(s))
end

# We want to solve the subproblems:
solutions = []
for s in 1:scenarios
    push!(solutions, Solve_Subproblem(subproblems[s]))
end

Then we obtein the error: OptimizeNotCalled()

I know I could join the creation and solving of subproblems in a single “for” loop, but the reality is that later in my code, it deals with modifying the subproblems by creating their augmented Lagrangian relaxations (very similar to the Progressive Hedging algorithm). So, in the iterations of my algorithm, I need to solve the subproblems multiple times, and it’s not optimal to have to create the subproblems for all scenarios in each iteration. I would like to be able to create them only once.

I hope I explained myself well. I would appreciate it if you could review the code. Thank you very much!

[benjamin.nunez]

Hello @mkitti, thank you very much for your response!

I will share with you the same code fragments that I shared with @odow:

Oscar (@odow) deleted the code because it was just a copy of the previous post

[odow]

first_stage_variables[1] = x

This stores the first stage variable x in the dictionary first_stage_variables with the key 1.

When you call Build_Subproblem with different values of s, this value gets overwritten. Thus, after Build_Subproblem, it contains the value of x from scenarios subproblem. When you solve the first subproblem, you query the value of the variable, but that problem hasn’t been optimized yet.

You need to think about a different data structure to use.

Here’s one approach:

using JuMP, Gurobi
scenarios = 3
probabiliity = [0.2, 0.5, 0.3]
demand = [[3.0, 3.0, 2.0], [5.0, 3.0, 2.0], [7.0, 3.0, 2.0]]
generators = 4 
periods = 3 
fixed_cost = [10, 7, 16, 6]
variable_cost = [40 24 4; 45 27 4.5; 32 19.2 3.2; 55 33 5.5]
minimum_power = 12
budget = 120
function build_subproblem(s::Int)
    subproblem = Model(Gurobi.Optimizer)
    set_silent(subproblem)
    @variable(subproblem, x[1:generators] >= 0)
    @variable(subproblem, y[1:generators, 1:periods] >= 0)
    @constraint(subproblem, sum(x) >= minimum_power)
    @constraint(subproblem, fixed_cost' * x <= budget)
    for i in 1:generators, j in 1:periods
        @constraint(subproblem, y[i, j] <= x[i])
    end
    for j in 1:periods
        @constraint(subproblem, sum(y[:, j]) >= demand[s][j])
    end
    @objective(subproblem, Min, fixed_cost' * x + sum(variable_cost .* y))
    return (model = subproblem, first_stage_variables = [x])
end

function solve_subproblem(data)
    optimize!(data.model)
    @assert is_solved_and_feasible(data.model)  # Important!!! Check solved okay
    first_stage_solutions = [
        JuMP.value.(variable) for variable in data.first_stage_variables
    ]
    return first_stage_solutions
end

subproblems = [build_subproblem(s) for s in 1:scenarios]
solutions = []
for s in 1:scenarios
    push!(solutions, solve_subproblem(subproblems[s]))
end

You will probably benefit from reading these tutorials:

• Design patterns for larger models · JuMP
• Example: two-stage newsvendor · SDDP.jl
• Example: deterministic to stochastic · SDDP.jl