stochastic variables make the SDDP model infeasible

[BaiYansong-bys]

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]

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]

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]

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

[BaiYansong-bys]

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

[odow]

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