odow / SDDP.jl

A JuMP extension for Stochastic Dual Dynamic Programming
https://sddp.dev
Other
293 stars 62 forks source link

stochastic variables make the SDDP model infeasible #640

Closed BaiYansong-bys closed 1 year ago

BaiYansong-bys commented 1 year ago

Hi, I have encountered such a problem, my model includes 20 stages, sce0 is vector{vector{Float64}} with 144 elements, Pro0 is vector{Float64} with 144 elements. pw_s is defined as stochastic variables. if I write this, the model can be solved:

if stage == 1
@constraint(sp, [i in 1:NW], pw[i]<=pw_l[i,stage]*(1)*uw[i])#pw_s[i]
else
SDDP.parameterize(sp,sce0,Pro0) do omga
     JuMP.fix.(pw_s,omga);
end
@constraint(sp, [i in 1:NW], pw[i]<=pw_l[i,stage]*(1)*uw[i])     
end

if I write this, the model cannot be solved:

if stage == 1
@constraint(sp, [i in 1:NW], pw[i]<=pw_l[i,stage]*(1)*uw[i])#pw_s[i]
else
SDDP.parameterize(sp,sce0,Pro0) do omga
     JuMP.fix.(pw_s,omga);
end
@constraint(sp, [i in 1:NW], pw[i]<=pw_l[i,stage]*(1+0*pw_s[i])*uw[i])     
end

The prompt error is as follows: ERROR: LoadError: Unable to retrieve solution from (20, 1). Termination status: OPTIMAL Primal status: FEASIBLE_POINT Dual status: NO_SOLUTION. A MathOptFormat file was written to subproblem.mof.json. See https://odow.github.io/SDDP.jl/latest/tutorial/06_warnings/#Numerical-stability-1 for more information. In theory, *pw_ s[i] 0==0**, why are the results different?

odow commented 1 year ago

Presumably, because your constraint is now a quadratic function, and the solver you are using cannot return the dual of a quadratic program.

I'd do something like

SDDP.parameterize(sp,sce0,Pro0) do omga
     for i in 1:NW
          # Note the -, needed because normalized constraint is pw[i] - uw[i] <= 0
          set_normalized_coefficient(c_stochastic, uw[i], -pw_l[i,stage]*(1+omga))
     end
end
@constraint(sp, c_stochastic[i in 1:NW], pw[i] <= uw[i]) 
BaiYansong-bys commented 1 year ago

Thanks a lot, Oscar, there still have some trouble. If I express the code as follows,

NW = 2
SDDP.parameterize(sp,[[0.2,0.3],[0,-0.1]],[0.5,0.5]) do omga
     for i in 1:NW
          # Note the -, needed because normalized constraint is pw[i] - uw[i] <= 0
          set_normalized_coefficient(c_stochastic, uw[i], -pw_l[i,stage]*(1+omga[i]))
     end
end
@constraint(sp, c_stochastic[i in 1:NW], pw[i] <= uw[i]) 

I find some errors, set_normalized_coefficient seems to have errors because the input cannot be "c_stochastic", the error is as follows: LoadError: MethodError: no method matching set_normalized_coefficient(::Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}, ::VariableRef, ::Float64) Closest candidates are: set_normalized_coefficient(::ConstraintRef{var"#s140", MathOptInterface.ConstraintIndex{F, S}, Shape} where {var"#s140"<:AbstractModel, Shape<:AbstractShape}, ::Any, ::Any) where {S, T, F<:Union{MathOptInterface.ScalarAffineFunction{T}, MathOptInterface.ScalarQuadraticFunction{T}}} at C:\Users\Lenovo.julia\packages\JuMP\klrjG\src\constraints.jl:596 I found this expression seems to work,

        SDDP.parameterize(sp, [[0.2,0.3],[0,-0.1]],[0.5,0.5]) do omga#sce0,Pro0
                #JuMP.fix.(pw_s, omga)
            #for i in 1:NW
                set_normalized_coefficient(c_stochastic1, uw[1], -pw_l[1,stage]*(1+omga[1]))
                set_normalized_coefficient(c_stochastic2, uw[2], -pw_l[2,stage]*(1+omga[2]))
                #
            #end
        end
        @constraint(sp, c_stochastic1, pw[1]<=uw[1])
        @constraint(sp, c_stochastic2, pw[2]<=uw[2])

Can the code above be more concise?

odow commented 1 year ago

I guess I meant set_normalized_coefficient(c_stochastic[i], I can't test the code without a reproducible example.

BaiYansong-bys commented 1 year ago

Thanks, Oscar, it can work with the above code.

odow commented 1 year ago

Closing as completed. Please re-open if you have any other questions.