Variable index allocation problem for silent variables

[ChrisFuelOR]

Hi,

currently I observe some issues with the variable allocation using GAMS.jl.

It seems that variables which are added to a JuMP model, when solving this model via GAMS.jl are only considered in the corresponding GAMS model if they do actually appear in one of the constraints (or the objective). While this may be intended to keep the GAMS models as small as possible, it seems to have some side effects on the index allocation. It seems that after solving the model, the variable values are not allocated to the original indices, but simply starting from one, without considering the "silent" variables, which have been omitted in GAMS.

Here is a minimal working example:

using GAMS
using JuMP

dual_vars = zeros(8)

approx_model = JuMP.Model(GAMS.Optimizer)
JuMP.set_optimizer_attribute(approx_model, "Solver", "Gurobi")

@variables approx_model begin
    θ
    x[1:length(dual_vars)]
end
JuMP.@objective(approx_model, JuMP.MOI.MAX_SENSE, θ)

JuMP.@constraint(
    approx_model,
    θ - 0.65 * x[3] + x[4] - x[5] + x[6] + x[7] <= 36
)

JuMP.optimize!(approx_model)
dual_vars .= value.(x)

At the last line, this yields "ERROR: LoadError: BoundsError: attempt to access 6-element Array{Float64,1} at index [7]".

If you look at the single outputs

println(JuMP.value(x[1]))
println(JuMP.value(x[2]))
println(JuMP.value(x[3]))
println(JuMP.value(x[4]))
println(JuMP.value(x[5]))
println(JuMP.value(x[6]))
println(JuMP.value(x[7]))
println(JuMP.value(x[8]))

you get

1.0
0.0
0.0
0.0
0.0
ERROR: LoadError: BoundsError: attempt to access 6-element Array{Float64,1} at index [7]

so it looks like the values of the variables x[3], x[4], x[5], x[6], x[7] appearing in the constraints are allocated to x[1] to x[5] now and for the silent variables no values are returned at all.

A bit more background where this issue occurs:

I try to apply a cutting-plane method to solve a Lagrangian dual problem (as for SDDiP in SDDP.jl), solving all subproblems using GAMS.jl. The dimension of the dual variables is determined right in the beginning. Let's say it's 8 as in the MWE, so 8 dual variables are added to the approx_model, which is then iteratively extended by cuts.

Now, if by solving the Lagrangian relaxation we obtain a subgradient with some 0 components in the first iteration, a cut is constructed, that contains only some of the dual variables. Solving approx_model considers only those variables then.

Updating dual_vars with the obtained dual solution yields the described error.

Many thanks in advance for having a look at this.

[odow]

Just chiming in here: I don't know about the GAMS issue, but please talk to me if you have development plans for SDDiP https://github.com/odow/SDDP.jl/issues/246!

[renkekuhlmann]

Thanks very much for reporting! I can confirm that the used GAMS savepoint option only includes the symbols that are present in the model, which results in the incorrect behaviour reported by you.

[renkekuhlmann]

This has been fixed with GAMS.jl v0.1.6. The output of the above program will be:

8-element Array{Float64,1}:
 NaN
 NaN
   1.0
   0.0
   0.0
   0.0
   0.0
 NaN

I thought about using missing instead of NaN, but that resulted in errors with MathOptInterface. missing would have been nice, but NaN works as well. Thanks again for reporting the issue!

[odow]

Isn't the correct answer 0, not NaN?

[renkekuhlmann]

Good question. I have thought about this for a while and I think that missing or NaN is best here as it indicates the value could have been set to anything since it wasn't part of the optimization. If I set it to something, e.g. 0, this is unclear to the user. But I am happy to hear your thoughts about it. Is there any convention in MathOptInterface to set the value to 0 in this case?

[odow]

I think the convention in most (every?) solver is that unspecified variables default to 0. That's what pure-Gurobi does, so it makes sense that is what GAMS-Gurobi should do:

using JuMP, Gurobi
approx_model = Model(Gurobi.Optimizer)
@variables(approx_model, begin
    θ <= 100
    x[1:8]
end)
@objective(approx_model, Max, θ)
@constraint(approx_model, θ - 0.65 * x[3] + x[4] - x[5] + x[6] + x[7] <= 36)
optimize!(approx_model)
julia> value.(x)
8-element Array{Float64,1}:
  0.0
  0.0
 98.46153846153845
  0.0
  0.0
  0.0
  0.0
  0.0

For SDDP, we take the solution to the Lagrangian dual, and use the optimal x as an input to a different LP. Getting NaN solutions is going to cause no-end of issues.

[renkekuhlmann]

I still doubt that 0 is "the correct answer", but

Getting NaN solutions is going to cause no-end of issues.

is a very good point. I will change the behaviour to return 0 with the next release. Thanks a lot for pointing this out!

[odow]

Any solution is correct. 0 is just the least worst option, and most often the one the user intended.

[ChrisFuelOR]

Thanks for your quick responses and for fixing the bug.