Handling the case when the subproblem is unbounded

[hideakiv]

Conditioned when the coupling variables are free, e.g., during the temporal decomposition

[anhphuong-ngo]

Hi @hideakiv and author, actually I have a similar problem. I develop a two-stage stochastic mixed-integer program, which is similar to the "sslp.jl" example of this repository, but the first stage of my model has 3 different sets of variables:

@variable(model, z[1:I], Bin)
@variable(model, y >= 0)
@variable(model, x[1:I]>=0)

The result is not the same as my expectation:

DD.Primal: Inf
DD.Dual: 4400

I tested the same model with other algorithms and direct solve (without using any decomposition methods), and they have quite similar results.

For DD algorithm, I checked the results shown on Julia terminal, at initial iterations the objective value is very close to the result of other algorithms. However, in the end, the last subproblems show no solution. When it stopped, DD.primal showed inf. I don't know what happened. That could be from coupling_variables (I raised it at question #52 ). It would be great if we can adjust any parameters or can fix it. Hope to hear from you. Thank you.

[hideakiv]

Hi @anhphuong-ngo, could you post the entire file so we can take a closer look?

[anhphuong-ngo]

Hi @hideakiv, here is mine:

using Pkg
Pkg.activate("SSLP_plus")
using DualDecomposition
using JuMP, Ipopt, GLPK, Gurobi
const DD = DualDecomposition

# Set and indices
I = 100
J = 20
NS = 20 
# Params
...

# This creates a Lagrange dual problem for each scenario s.
function create_scenario_model(s::Int64)
    # m = Model(GLPK.Optimizer)
    m = Model(Gurobi.Optimizer)

    @variable(m, 0 <= xe[1:I,1:J] <= 100000000)            # 2nd var. The reason I set 0 <= xen <= 10000.., because when I set lower_bound = 0, it ran into error
    @variable(m, 0 <= xc[1:I,] <= 100000000)               # 2nd var
    @variable(m, 0 <= ye[1:J] <= 100000000)                # 1st var
    @variable(m, 0 <= y <= 100000000)                      # 1st var
    @variable(m, 0 <= z[1:J] <= 1, Bin)                    # 1st var 

...

    return m
end

# Create DualDecomposition instance.
algo = DD.LagrangeDual()

# Add Lagrange dual problem for each scenario s.
models = Dict{Int,JuMP.Model}(s => create_scenario_model(s) for s in 1:NS)
for s in 1:NS
    DD.add_block_model!(algo, s, models[s])
end

coupling_variables = Vector{DD.CouplingVariableRef}()
for s in 1:NS
    model = models[s]
    zref = model[:z]
    for j in 1:J
        push!(coupling_variables, DD.CouplingVariableRef(s, j, zref[j]))
    end

    ynref = model[:ye]
    for j in 1:J
        push!(coupling_variables, DD.CouplingVariableRef(s, j, ynref[j]))
    end

    push!(coupling_variables, DD.CouplingVariableRef(s, "y", model[:y]))
end

# Set nonanticipativity variables as an array of symbols.
DD.set_coupling_variables!(algo, coupling_variables)

# Lagrange master method
# LM = DD.BundleMaster(BM.ProximalMethod, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
LM = DD.BundleMaster(BM.ProximalMethod, optimizer_with_attributes(Gurobi.Optimizer))

# add heuristic
DD.add!(DD.AllBlockHeuristic, algo)
DD.add!(DD.RoundingHeuristic, algo)

# Solve the problem with the solver; this solver is for the underlying bundle method.
DD.run!(algo, LM)

@show DD.primal_objective_value(algo)
@show DD.dual_objective_value(algo)
@show DD.primal_solution(algo)
@show DD.dual_solution(algo)

It would be great if you can have a look. Thanks.

[hideakiv]

I am getting infeasible status for the subproblems... did you encounter the same issue?

[hideakiv]

It appears that you are assigning the same coupling_id for z[j] and ye[j]. Could you try the following? push!(coupling_variables, DD.CouplingVariableRef(s, "z_$j", zref[j])) push!(coupling_variables, DD.CouplingVariableRef(s, "yn_$j", ynref[j]))

Although I doubt it resolves the subproblem issue...

[anhphuong-ngo]

I am getting infeasible status for the subproblems... did you encounter the same issue?

Thanks for your reply. Here is the result I recieved:

User-callback calls 31, time in user-callback 0.00 sec
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 4 rows, 821 columns and 61 nonzeros
Model fingerprint: 0xc803b2d3
Model has 820 quadratic objective terms
Coefficient statistics:
  Matrix range     [1e+00, 2e+01]
  Objective range  [1e+00, 1e+00]
  QObjective range [1e-02, 1e-02]
  Bounds range     [0e+00, 0e+00]
  RHS range        [0e+00, 0e+00]
Presolve removed 1 rows and 761 columns
Presolve time: 0.00s
Presolved: 3 rows, 60 columns, 60 nonzeros
Presolved model has 60 quadratic objective terms
Ordering time: 0.00s

Barrier statistics:
 Free vars  : 60
 AA' NZ     : 0.000e+00
 Factor NZ  : 3.000e+00
 Factor Ops : 3.000e+00 (less than 1 second per iteration)
 Threads    : 1

Barrier solved model in 0 iterations and 0.00 seconds (0.00 work units)
Optimal objective 0.00000000e+00

User-callback calls 43, time in user-callback 0.00 sec
Iter    1: ncols   821, nrows     4, fx0 -2.000000e+01, fy -2.000000e+01, m -2.000000e+01, v 0.000000e+00, u 1.000000e-02, i +0, master time    0.0s, eval time    0.7s, time    0.8s
TERMINATION: Optimal: v = 0.0
DD.primal_objective_value(algo) = Inf
DD.dual_objective_value(algo) = 20.0
DD.primal_solution(algo) = Dict{Int64, Float64}()

[anhphuong-ngo]

It appears that you are assigning the same coupling_id for z[j] and ye[j]. Could you try the following? push!(coupling_variables, DD.CouplingVariableRef(s, "z_$j", zref[j])) push!(coupling_variables, DD.CouplingVariableRef(s, "yn_$j", ynref[j]))

Although I doubt it resolves the subproblem issue...

Thanks for your reply @hideakiv . I tried with your suggestions, but it seems subproblem could fina a solution, but the it does not convert:

...
Optimal solution found (tolerance 1.00e-04)
Best objective 1.000000000000e+00, best bound 1.000000000000e+00, gap 0.0000%

User-callback calls 25, time in user-callback 0.00 sec
ERROR: MethodError: Cannot `convert` an object of type String to an object of type Int64

Closest candidates are:
  convert(::Type{T}, ::Base.TwicePrecision) where T<:Number
   @ Base twiceprecision.jl:273
  convert(::Type{T}, ::AbstractChar) where T<:Number
   @ Base char.jl:185
  convert(::Type{T}, ::CartesianIndex{1}) where T<:Number
   @ Base multidimensional.jl:127
  ...

Stacktrace:
  [1] setindex!(h::Dict{Int64, Float64}, v0::Float64, key0::String)
    @ Base .\dict.jl:361
  [2] Dict{Int64, Float64}(kv::Dict{Any, Any})
    @ Base .\dict.jl:84
  [3] convert
    @ .\abstractdict.jl:568 [inlined]
  [4] setproperty!
    @ .\Base.jl:38 [inlined]
  [5] run!(#unused#::Type{DualDecomposition.AllBlockHeuristic}, LD::DualDecomposition.LagrangeDual, val::Dict{Int64, SparseArrays.SparseVector{Float64}}, ub::Dict{Int64, SparseArrays.SparseVector{Float64}}, lb::Dict{Int64, SparseArrays.SparseVector{Float64}})
    @ DualDecomposition C:\Users\ango1\.julia\packages\DualDecomposition\rqrRm\src\PrimalHeuristics.jl:120
  [6] (::DualDecomposition.var"#solveLagrangeDual#9"{DualDecomposition.LagrangeDual, Int64})(λ::Vector{Float64})
    @ DualDecomposition C:\Users\ango1\.julia\packages\DualDecomposition\rqrRm\src\LagrangeDual.jl:161
  [7] evaluate_functions!(method::BundleMethod.ProximalMethod)
    @ BundleMethod C:\Users\ango1\.julia\packages\BundleMethod\H864L\src\Proximal.jl:160
  [8] add_initial_bundles!
    @ C:\Users\ango1\.julia\packages\BundleMethod\H864L\src\Abstract.jl:67 [inlined]
  [9] run!(method::BundleMethod.ProximalMethod)
    @ BundleMethod C:\Users\ango1\.julia\packages\BundleMethod\H864L\src\Abstract.jl:52
 [10] run!(method::DualDecomposition.BundleMaster)
    @ DualDecomposition C:\Users\ango1\.julia\packages\DualDecomposition\rqrRm\src\LagrangeMaster\BundleMethod.jl:52
 [11] run!(LD::DualDecomposition.LagrangeDual, LM::DualDecomposition.BundleMaster, initial_λ::Nothing)
    @ DualDecomposition C:\Users\ango1\.julia\packages\DualDecomposition\rqrRm\src\LagrangeDual.jl:200
 [12] run!(LD::DualDecomposition.LagrangeDual, LM::DualDecomposition.BundleMaster)
    @ DualDecomposition C:\Users\ango1\.julia\packages\DualDecomposition\rqrRm\src\LagrangeDual.jl:69
 [13] top-level scope
    @ c:\Users\ango1\MyDSPOpt\Testing_sslp_plus.jl:94

[hideakiv]

Ah then would you mind changing it to push!(coupling_variables, DD.CouplingVariableRef(s, j, zref[j])) push!(coupling_variables, DD.CouplingVariableRef(s, J+j, ynref[j])) push!(coupling_variables, DD.CouplingVariableRef(s, 0, model[:y]))

[anhphuong-ngo]

@hideakiv Thanks for your help. It works well. Just one more question. For this one push!(coupling_variables, DD.CouplingVariableRef(s, 0, model[:y]))

I tried with push!(coupling_variables, DD.CouplingVariableRef(s, 0, model[:y])) and push!(coupling_variables, DD.CouplingVariableRef(s, 2*J+1, model[:y]))

Both are working, but according to your experience, which one is better, indexing starting from 0 or 1? Thanks.

[hideakiv]

I would suggest a readable coupling ID, like strings. This should be fixed in #53.