Conic subproblem numerical issue

[haoxiangyang89]

I am running a conic subproblem for SOCP relaxation of the OPF problem. For the first a few iteration of SDDP, it works fine. However, for the test case I had, at about iteration 25, I saw Gurobi error 10005. The error message is displayed below:

SDDP.train(model; stopping_rules = [SDDP.BoundStalling(10, 1e-4)], duality_handler = SDDP.ContinuousConicDuality())

                      SDDP.jl (c) Oscar Dowson, 2017-21

Problem
  Nodes           : 8
  State variables : 19
  Scenarios       : 6.56100e+03
  Existing cuts   : false
  Subproblem structure                                              : (min, max)
    Variables                                                       : (202, 202)
    GenericAffExpr{Float64,VariableRef} in MOI.GreaterThan{Float64} : (39, 39)
    Array{VariableRef,1} in MOI.SecondOrderCone                     : (31, 31)
    VariableRef in MOI.LessThan{Float64}                            : (27, 28)
    VariableRef in MOI.GreaterThan{Float64}                         : (80, 80)
    GenericAffExpr{Float64,VariableRef} in MOI.LessThan{Float64}    : (67, 67)
    GenericAffExpr{Float64,VariableRef} in MOI.EqualTo{Float64}     : (105, 112)
Options
  Solver          : serial mode
  Risk measure    : SDDP.Expectation()
  Sampling scheme : SDDP.InSampleMonteCarlo

Numerical stability report
  Non-zero Matrix range     [8e-04, 2e+00]
  Non-zero Objective range  [0e+00, 0e+00]
  Non-zero Bounds range     [9e-01, 1e+00]
  Non-zero RHS range        [1e-02, 9e+03]
No problems detected

 Iteration    Simulation       Bound         Time (s)    Proc. ID   # Solves
        1    1.844339e+04   1.841453e+04   2.881200e-01          1         32
        2    1.655765e+04   1.841870e+04   6.299469e-01          1         64
        3    1.655285e+04   1.841871e+04   8.955100e-01          1         96
        4    1.655186e+04   1.841872e+04   1.156718e+00          1        128
        5    1.655186e+04   1.841872e+04   1.407547e+00          1        160
        6    1.655186e+04   1.841872e+04   1.667418e+00          1        192
        7    1.655186e+04   1.841872e+04   1.921739e+00          1        224
        8    1.655187e+04   1.841872e+04   2.178782e+00          1        256
        9    1.655187e+04   1.841870e+04   2.434639e+00          1        288
       10    1.655185e+04   1.841870e+04   2.742453e+00          1        320
       11    1.655185e+04   1.841871e+04   3.012145e+00          1        352
       12    1.655186e+04   1.841871e+04   3.270867e+00          1        384
       13    1.655186e+04   1.841871e+04   3.536912e+00          1        416
       14    1.655186e+04   1.841871e+04   3.817132e+00          1        448
       15    1.655186e+04   1.841871e+04   4.093773e+00          1        480
       16    1.655186e+04   1.841872e+04   4.374855e+00          1        512
       17    1.655186e+04   1.841870e+04   4.644258e+00          1        544
       18    1.655185e+04   1.841870e+04   4.950060e+00          1        576
       19    1.655185e+04   1.841870e+04   5.247546e+00          1        608
       20    1.655185e+04   1.841872e+04   5.561274e+00          1        640
       21    1.655186e+04   1.841870e+04   5.840731e+00          1        672
       22    1.655185e+04   1.841872e+04   6.151065e+00          1        704
       23    1.655186e+04   1.841870e+04   6.432530e+00          1        736
       24    1.655185e+04   1.841870e+04   6.743254e+00          1        768
       25    1.655185e+04   1.841871e+04   7.045040e+00          1        800
ERROR: Gurobi Error 10005: Unable to retrieve attribute 'RC'
Stacktrace:
 [1] _check_ret at /Users/haoxiangyang/.julia/packages/Gurobi/BAtIN/src/MOI_wrapper/MOI_wrapper.jl:326 [inlined]
 [2] get(::Gurobi.Optimizer, ::MathOptInterface.ConstraintDual, ::MathOptInterface.ConstraintIndex{MathOptInterface.SingleVariable,MathOptInterface.EqualTo{Float64}}) at /Users/haoxiangyang/.julia/packages/Gurobi/BAtIN/src/MOI_wrapper/MOI_wrapper.jl:3023
 [3] get(::MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, ::MathOptInterface.ConstraintDual, ::MathOptInterface.ConstraintIndex{MathOptInterface.SingleVariable,MathOptInterface.EqualTo{Float64}}) at /Users/haoxiangyang/.julia/packages/MathOptInterface/YDdD3/src/Bridges/bridge_optimizer.jl:1208
 [4] get(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64,MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}}, ::MathOptInterface.ConstraintDual, ::MathOptInterface.ConstraintIndex{MathOptInterface.SingleVariable,MathOptInterface.EqualTo{Float64}}) at /Users/haoxiangyang/.julia/packages/MathOptInterface/YDdD3/src/Utilities/cachingoptimizer.jl:757
 [5] _moi_get_result(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64,MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}}, ::MathOptInterface.ConstraintDual, ::Vararg{Any,N} where N) at /Users/haoxiangyang/.julia/packages/JuMP/klrjG/src/JuMP.jl:1199
 [6] get(::Model, ::MathOptInterface.ConstraintDual, ::ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.SingleVariable,MathOptInterface.EqualTo{Float64}},ScalarShape}) at /Users/haoxiangyang/.julia/packages/JuMP/klrjG/src/JuMP.jl:1244
 [7] _constraint_dual(::ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.SingleVariable,MathOptInterface.EqualTo{Float64}},ScalarShape}, ::Int64) at /Users/haoxiangyang/.julia/packages/JuMP/klrjG/src/constraints.jl:865
 [8] #186 at /Users/haoxiangyang/.julia/packages/JuMP/klrjG/src/constraints.jl:851 [inlined]
 [9] iterate at ./generator.jl:47 [inlined]
 [10] Dict{Symbol,Float64}(::Base.Generator{Dict{Symbol,SDDP.State{VariableRef}},SDDP.var"#186#187"{Int64}}) at ./dict.jl:102
 [11] #solve_subproblem#61(::SDDP.ContinuousConicDuality, ::typeof(SDDP.solve_subproblem), ::SDDP.PolicyGraph{Int64}, ::SDDP.Node{Int64}, ::Dict{Symbol,Float64}, ::Int64, ::Array{Tuple{Int64,Any},1}) at /Users/haoxiangyang/.julia/packages/SDDP/LvBi8/src/plugins/duality_handlers.jl:124
 [12] (::SDDP.var"#kw##solve_subproblem")(::NamedTuple{(:duality_handler,),Tuple{SDDP.ContinuousConicDuality}}, ::typeof(SDDP.solve_subproblem), ::SDDP.PolicyGraph{Int64}, ::SDDP.Node{Int64}, ::Dict{Symbol,Float64}, ::Int64, ::Array{Tuple{Int64,Any},1}) at ./none:0
 [13] macro expansion at /Users/haoxiangyang/.julia/packages/SDDP/LvBi8/src/algorithm.jl:663 [inlined]
 [14] macro expansion at /Users/haoxiangyang/.julia/packages/TimerOutputs/SSeq1/src/TimerOutput.jl:236 [inlined]
 [15] solve_all_children(::SDDP.PolicyGraph{Int64}, ::SDDP.Node{Int64}, ::SDDP.BackwardPassItems{Int64,SDDP.Noise}, ::Float64, ::Nothing, ::Nothing, ::Dict{Symbol,Float64}, ::SDDP.CompleteSampler, ::Array{Tuple{Int64,Any},1}, ::SDDP.ContinuousConicDuality) at /Users/haoxiangyang/.julia/packages/SDDP/LvBi8/src/algorithm.jl:662
 [16] backward_pass(::SDDP.PolicyGraph{Int64}, ::SDDP.Options{Int64}, ::Array{Tuple{Int64,Any},1}, ::Array{Dict{Symbol,Float64},1}, ::Array{Tuple{},1}, ::Array{Tuple{Int64,Dict{Int64,Float64}},1}) at /Users/haoxiangyang/.julia/packages/SDDP/LvBi8/src/algorithm.jl:532
 [17] macro expansion at /Users/haoxiangyang/.julia/packages/SDDP/LvBi8/src/algorithm.jl:782 [inlined]
 [18] macro expansion at /Users/haoxiangyang/.julia/packages/TimerOutputs/SSeq1/src/TimerOutput.jl:236 [inlined]
 [19] iteration(::SDDP.PolicyGraph{Int64}, ::SDDP.Options{Int64}) at /Users/haoxiangyang/.julia/packages/SDDP/LvBi8/src/algorithm.jl:781
 [20] master_loop(::SDDP.Serial, ::SDDP.PolicyGraph{Int64}, ::SDDP.Options{Int64}) at /Users/haoxiangyang/.julia/packages/SDDP/LvBi8/src/plugins/parallel_schemes.jl:32
 [21] #train#67(::Nothing, ::Nothing, ::Int64, ::String, ::Int64, ::Bool, ::Array{SDDP.BoundStalling,1}, ::SDDP.Expectation, ::SDDP.InSampleMonteCarlo, ::SDDP.CutType, ::Float64, ::Bool, ::Int64, ::SDDP.CompleteSampler, ::Bool, ::SDDP.Serial, ::SDDP.DefaultForwardPass, ::Nothing, ::Bool, ::SDDP.ContinuousConicDuality, ::SDDP.var"#71#75", ::typeof(SDDP.train), ::SDDP.PolicyGraph{Int64}) at /Users/haoxiangyang/.julia/packages/SDDP/LvBi8/src/algorithm.jl:1036
 [22] (::SDDP.var"#kw##train")(::NamedTuple{(:stopping_rules, :duality_handler),Tuple{Array{SDDP.BoundStalling,1},SDDP.ContinuousConicDuality}}, ::typeof(SDDP.train), ::SDDP.PolicyGraph{Int64}) at ./none:0
 [23] top-level scope at REPL[130]:1

I have turned on the QCPDual option for Gurobi so it can generate conic dual. I am just wondering if this is a general issue or is there any way to improve the stability.

Also I checked the solution output by the simulation results. It seems for the constraints are changed and all scenarios need to follow those constraints. What I would like to achieve is to have scenario specific constraints. According to different contingencies, different sets of constraints are applied.

The test can be found here: https://github.com/haoxiangyang89/disruptionN-1/blob/master/src/convex/ms_sddp.jl. @odow should have access to this repo. Thanks!

[odow]

Ah nice! I've been trying to find a test-case for this: https://github.com/jump-dev/Gurobi.jl/issues/415.

Can you run it to get the Gurobi log?

You can save the file with SDDP.write_to_file(model, "gurobi_failure.mof.json.gz"; test_scenarios = 0)

I assume what happens is that the solution is at the point of the cone where the dual is not defined.

What happens if you add sum(x) >= 0.0001 for your x in SecondOrderCone() constraints?

[odow]

It seems for the constraints are changed and all scenarios need to follow those constraints. What I would like to achieve is to have scenario specific constraints. According to different contingencies, different sets of constraints are applied.

Your parameterize code is incorrect:

        SDDP.parameterize(subproblem, support, nominal_probability) do ω
            if ω isa Int
                # Change generator bounds
                for g in fData.genIDList
                    if g == ω
                        JuMP.set_normalized_rhs(spUb[ω], 0.0)
                        JuMP.set_normalized_rhs(sqUb[ω], 0.0)
                        JuMP.set_normalized_rhs(spLb[ω], 0.0)
                        JuMP.set_normalized_rhs(sqLb[ω], 0.0)
                        JuMP.set_normalized_rhs(genRamp_up[ω], 1000*fData.RU[g])
                        JuMP.set_normalized_rhs(genRamp_down[ω], 1000*fData.RD[g])
                    end
                end
            else
                # Change line bounds
                for l in fData.brList
                    if ((l[1] == ω[1])&&(l[2] == ω[2]))||((l[1] == ω[2])&&(l[2] == ω[1]))
                        JuMP.set_normalized_rhs(linDistFlow_ub[l], -1000*fData.rateA[l]);
                        JuMP.set_normalized_rhs(linDistFlow_lb[l], 1000*fData.rateA[l]);
                        JuMP.set_normalized_rhs(thermal[l], 0.0);
                    end
                end
            end
        end

It's not sufficient to just change the current values, you need to set the values for all elements that change in any of the scenarios (i.e., we don't set values back to their defaults).

[haoxiangyang89]

It seems for the constraints are changed and all scenarios need to follow those constraints. What I would like to achieve is to have scenario specific constraints. According to different contingencies, different sets of constraints are applied.

Your parameterize code is incorrect:

        SDDP.parameterize(subproblem, support, nominal_probability) do ω
            if ω isa Int
                # Change generator bounds
                for g in fData.genIDList
                    if g == ω
                        JuMP.set_normalized_rhs(spUb[ω], 0.0)
                        JuMP.set_normalized_rhs(sqUb[ω], 0.0)
                        JuMP.set_normalized_rhs(spLb[ω], 0.0)
                        JuMP.set_normalized_rhs(sqLb[ω], 0.0)
                        JuMP.set_normalized_rhs(genRamp_up[ω], 1000*fData.RU[g])
                        JuMP.set_normalized_rhs(genRamp_down[ω], 1000*fData.RD[g])
                    end
                end
            else
                # Change line bounds
                for l in fData.brList
                    if ((l[1] == ω[1])&&(l[2] == ω[2]))||((l[1] == ω[2])&&(l[2] == ω[1]))
                        JuMP.set_normalized_rhs(linDistFlow_ub[l], -1000*fData.rateA[l]);
                        JuMP.set_normalized_rhs(linDistFlow_lb[l], 1000*fData.rateA[l]);
                        JuMP.set_normalized_rhs(thermal[l], 0.0);
                    end
                end
            end
        end

It's not sufficient to just change the current values, you need to set the values for all elements that change in any of the scenarios (i.e., we don't set values back to their defaults).

Thanks! I fixed the rhs part. The Gurobi error disappears when we add a strictly positive constraints to the second order conic constraint, sum(|x_i|) >= 0.00001. I am not sure if I am running it correctly, but after I train the SDDP and saw the error of Gurobi 10005, I tried to do SDDP.write_to_file(model, "gurobi_failure.mof.json.gz"; test_scenarios = 0)

I got this error: ERROR: StochOptFormat does not support writing after a call to 'SDDP.train'.

But when I only set up the model and run the write_to_file, it gives me this error: ERROR: MethodError: Cannot 'convert' an object of type GenericQuadExpr{Float64,VariableRef} to an object of type GenericAffExpr{Float64,VariableRef}

I have commited the change to my repo. The added constraints are in Line 128-142. Commenting them out will give the error of Gurobi 10005.

[haoxiangyang89]

I thought about this a little bit more. When we generate the cut, sddp.jl uses the dual value obtained from the solver. Is there an option to derive the dual problem, solve the dual, and then directly use the dual variable's value obtained from this dual problem? My past experience is that when we directly derive the dual (I performed this step by hand), at least we can get a dual value and it is usually much more numerically stable than relying on the dual output from the solver.

[odow]

it gives me this error:

Yeah you need to call write_to_file before train. I'll take a look at the quadratic error. Do you have the full stack trace?

The Gurobi error disappears when we add a strictly positive constraints to the second order conic constraint

The dual is not defined at the point of the cone, because you have a 0 >= sqrt(0) issue.

Is there an option to derive the dual problem, solve the dual, and then directly use the dual variable's value obtained from this dual problem

No. When we encounter issues like this we usually throw away the basis matrix and start a fresh solve. The issue is that Gurobi.jl lies to us and says it has a dual solution when it actually doesn't.

[odow]

@haoxiangyang89 how do I reproduce this? I tried running the ms_sddp.jl file, and I got ci not defined.

[odow]

I tried with ci=1 and the constraints you added removed, but it runs okay. Also: what is going on with your bound? It's not monotonic. Try with cut_deletion_minimum = 10_000.

[odow]

Closing in favor of https://github.com/jump-dev/Gurobi.jl/issues/415. This doesn't seem to be something we can fix in SDDP.jl.