martinbiel / StochasticPrograms.jl

Julia package for formulating and analyzing stochastic recourse models.
MIT License
75 stars 25 forks source link

Can I use StochasticPrograms.jl with Conic Constraints? #39

Open jmcastro2109 opened 2 years ago

jmcastro2109 commented 2 years ago

Hi,

My second stage problem has a linear objective and the constraint is a power cone. Currently I am solving this second stage problem with Mosek. Can I use Stochastic Programs to solve the first stage?

edit:

Looks like it does not work :(

Solver name: Mosek
ERROR: LoadError: Unable to transform a quadratic constraint into a second-order cone
constraint because the quadratic constraint is not strongly convex.

Convex constraints that are not strongly convex (i.e., the matrix is
positive semidefinite but not positive definite) are not supported
yet.

Note that a quadratic equality constraint is non-convex.
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:33
  [2] bridge_constraint(#unused#::Type{MathOptInterface.Bridges.Constraint.QuadtoSOCBridge{Float64}}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, func::MathOptInterface.ScalarQuadraticFunction{Float64}, set::MathOptInterface.GreaterThan{Float64})
    @ MathOptInterface.Bridges.Constraint C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\Constraint\quad_to_soc.jl:82
  [3] add_bridged_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, BridgeType::Type, f::MathOptInterface.ScalarQuadraticFunction{Float64}, s::MathOptInterface.GreaterThan{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:1404
  [4] add_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, f::MathOptInterface.ScalarQuadraticFunction{Float64}, s::MathOptInterface.GreaterThan{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:1483
  [5] normalize_and_add_constraint(model::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, func::MathOptInterface.ScalarQuadraticFunction{Float64}, set::MathOptInterface.GreaterThan{Float64}; allow_modify_function::Bool)
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\constraints.jl:19
  [6] normalize_and_add_constraint(model::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, func::MathOptInterface.ScalarQuadraticFunction{Float64}, set::MathOptInterface.GreaterThan{Float64})
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\constraints.jl:19
  [7] bridge_objective(#unused#::Type{MathOptInterface.Bridges.Objective.SlackBridge{Float64, MathOptInterface.ScalarQuadraticFunction{Float64}, MathOptInterface.ScalarQuadraticFunction{Float64}}}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, func::MathOptInterface.ScalarQuadraticFunction{Float64})
    @ MathOptInterface.Bridges.Objective C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\Objective\slack.jl:43
  [8] _bridge_objective(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, BridgeType::Type, func::MathOptInterface.ScalarQuadraticFunction{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:957
  [9] set(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, attr::MathOptInterface.ObjectiveFunction{MathOptInterface.ScalarQuadraticFunction{Float64}}, func::MathOptInterface.ScalarQuadraticFunction{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:995
 [10] bridge_objective
    @ C:\Users\juan\.julia\packages\StochasticPrograms\Jl6sf\src\types\decisions\bridges\objectives\quadratic.jl:13 [inlined]
 [11] _bridge_objective(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, BridgeType::Type, func::QuadraticDecisionFunction{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:957
 [12] set(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, attr::MathOptInterface.ObjectiveFunction{QuadraticDecisionFunction{Float64}}, func::QuadraticDecisionFunction{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:995
 [13] _pass_attributes(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}, idxmap::MathOptInterface.Utilities.IndexMap, attrs::Vector{MathOptInterface.AbstractModelAttribute}, supports_args::Tuple{}, get_args::Tuple{}, set_args::Tuple{}, pass_attr!::Function)
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:300
 [14] pass_attributes
    @ C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:206 [inlined]
 [15] pass_attributes
    @ C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:202 [inlined]
 [16] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}, copy_names::Bool, filter_constraints::Nothing)
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:711
 [17] #automatic_copy_to#127
    @ C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:24 [inlined]
 [18] #copy_to#4
    @ C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:401 [inlined]
 [19] attach_optimizer(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}})
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\cachingoptimizer.jl:185
 [20] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}})
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\cachingoptimizer.jl:248
 [21] optimize!(structure::DeterministicEquivalent{2, 1, Tuple{Vector{Scenario{NamedTuple{(:z₁, :z₂), Tuple{Float64, Float64}}}}}}, optimizer::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}}, x₀::Vector{Float64})
    @ StochasticPrograms C:\Users\juan\.julia\packages\StochasticPrograms\Jl6sf\src\methods\deterministic_equivalent\optimization.jl:13
 [22] optimize!(stochasticprogram::StochasticProgram{2, Tuple{StochasticPrograms.Stage{NamedTuple{(), Tuple{}}}, StochasticPrograms.Stage{NamedTuple{(), Tuple{}}}}, DeterministicEquivalent{2, 1, Tuple{Vector{Scenario{NamedTuple{(:z₁, :z₂), Tuple{Float64, Float64}}}}}}}; crash::StochasticPrograms.Crash.None, cache::Bool, kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ StochasticPrograms C:\Users\juan\.julia\packages\StochasticPrograms\Jl6sf\src\methods\api.jl:209
 [23] optimize!(stochasticprogram::StochasticProgram{2, Tuple{StochasticPrograms.Stage{NamedTuple{(), Tuple{}}}, StochasticPrograms.Stage{NamedTuple{(), Tuple{}}}}, DeterministicEquivalent{2, 1, Tuple{Vector{Scenario{NamedTuple{(:z₁, :z₂), Tuple{Float64, Float64}}}}}}})
    @ StochasticPrograms C:\Users\juan\.julia\packages\StochasticPrograms\Jl6sf\src\methods\api.jl:201
 [24] top-level scope
    @ c:\Users\juan\Dropbox\Projects\JMP\others\MWE_2.jl:196
in expression starting at c:\Users\juan\Dropbox\Projects\JMP\others\MWE_2.jl:196

caused by: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
  [1] cholesky!(F::SuiteSparse.CHOLMOD.Factor{Float64}, A::SuiteSparse.CHOLMOD.Sparse{Float64}; shift::Float64, check::Bool)
    @ SuiteSparse.CHOLMOD C:\Users\juan\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\SuiteSparse\src\cholmod.jl:1150
  [2] #cholesky#8
    @ C:\Users\juan\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\SuiteSparse\src\cholmod.jl:1185 [inlined]
  [3] cholesky
    @ C:\Users\juan\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\SuiteSparse\src\cholmod.jl:1182 [inlined]
  [4] #cholesky#9
    @ C:\Users\juan\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\SuiteSparse\src\cholmod.jl:1297 [inlined]
  [5] cholesky
    @ C:\Users\juan\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\SuiteSparse\src\cholmod.jl:1297 [inlined]
  [6] bridge_constraint(#unused#::Type{MathOptInterface.Bridges.Constraint.QuadtoSOCBridge{Float64}}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, func::MathOptInterface.ScalarQuadraticFunction{Float64}, set::MathOptInterface.GreaterThan{Float64})
    @ MathOptInterface.Bridges.Constraint C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\Constraint\quad_to_soc.jl:78
  [7] add_bridged_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, BridgeType::Type, f::MathOptInterface.ScalarQuadraticFunction{Float64}, s::MathOptInterface.GreaterThan{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:1404
  [8] add_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, f::MathOptInterface.ScalarQuadraticFunction{Float64}, s::MathOptInterface.GreaterThan{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:1483
  [9] normalize_and_add_constraint(model::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, func::MathOptInterface.ScalarQuadraticFunction{Float64}, set::MathOptInterface.GreaterThan{Float64}; allow_modify_function::Bool)
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\constraints.jl:19
 [10] normalize_and_add_constraint(model::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, func::MathOptInterface.ScalarQuadraticFunction{Float64}, set::MathOptInterface.GreaterThan{Float64})
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\constraints.jl:19
 [11] bridge_objective(#unused#::Type{MathOptInterface.Bridges.Objective.SlackBridge{Float64, MathOptInterface.ScalarQuadraticFunction{Float64}, MathOptInterface.ScalarQuadraticFunction{Float64}}}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, func::MathOptInterface.ScalarQuadraticFunction{Float64})
    @ MathOptInterface.Bridges.Objective C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\Objective\slack.jl:43
 [12] _bridge_objective(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, BridgeType::Type, func::MathOptInterface.ScalarQuadraticFunction{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:957
 [13] set(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, attr::MathOptInterface.ObjectiveFunction{MathOptInterface.ScalarQuadraticFunction{Float64}}, func::MathOptInterface.ScalarQuadraticFunction{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:995
 [14] bridge_objective
    @ C:\Users\juan\.julia\packages\StochasticPrograms\Jl6sf\src\types\decisions\bridges\objectives\quadratic.jl:13 [inlined]
 [15] _bridge_objective(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, BridgeType::Type, func::QuadraticDecisionFunction{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:957
 [16] set(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, attr::MathOptInterface.ObjectiveFunction{QuadraticDecisionFunction{Float64}}, func::QuadraticDecisionFunction{Float64})
    @ MathOptInterface.Bridges C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:995
 [17] _pass_attributes(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}, idxmap::MathOptInterface.Utilities.IndexMap, attrs::Vector{MathOptInterface.AbstractModelAttribute}, supports_args::Tuple{}, get_args::Tuple{}, set_args::Tuple{}, pass_attr!::Function)
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:300
 [18] pass_attributes
    @ C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:206 [inlined]
 [19] pass_attributes
    @ C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:202 [inlined]
 [20] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{MosekModel}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}, copy_names::Bool, filter_constraints::Nothing)
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:711
 [21] #automatic_copy_to#127
    @ C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\copy.jl:24 [inlined]
 [22] #copy_to#4
    @ C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Bridges\bridge_optimizer.jl:401 [inlined]
 [23] attach_optimizer(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}})
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\cachingoptimizer.jl:185
 [24] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}})
    @ MathOptInterface.Utilities C:\Users\juan\.julia\packages\MathOptInterface\YDdD3\src\Utilities\cachingoptimizer.jl:248
 [25] optimize!(structure::DeterministicEquivalent{2, 1, Tuple{Vector{Scenario{NamedTuple{(:z₁, :z₂), Tuple{Float64, Float64}}}}}}, optimizer::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}}, x₀::Vector{Float64})
    @ StochasticPrograms C:\Users\juan\.julia\packages\StochasticPrograms\Jl6sf\src\methods\deterministic_equivalent\optimization.jl:13
 [26] optimize!(stochasticprogram::StochasticProgram{2, Tuple{StochasticPrograms.Stage{NamedTuple{(), Tuple{}}}, StochasticPrograms.Stage{NamedTuple{(), Tuple{}}}}, DeterministicEquivalent{2, 1, Tuple{Vector{Scenario{NamedTuple{(:z₁, :z₂), Tuple{Float64, Float64}}}}}}}; crash::StochasticPrograms.Crash.None, cache::Bool, kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ StochasticPrograms C:\Users\juan\.julia\packages\StochasticPrograms\Jl6sf\src\methods\api.jl:209
 [27] optimize!(stochasticprogram::StochasticProgram{2, Tuple{StochasticPrograms.Stage{NamedTuple{(), Tuple{}}}, StochasticPrograms.Stage{NamedTuple{(), Tuple{}}}}, DeterministicEquivalent{2, 1, Tuple{Vector{Scenario{NamedTuple{(:z₁, :z₂), Tuple{Float64, Float64}}}}}}})
    @ StochasticPrograms C:\Users\juan\.julia\packages\StochasticPrograms\Jl6sf\src\methods\api.jl:201
 [28] top-level scope
    @ c:\Users\juan\Dropbox\Projects\JMP\others\MWE_2.jl:196
martinbiel commented 2 years ago

From the error text I would not assume that this is a StochasticPrograms related error. Does it work if you formulate a deterministic version of the problem in JuMP directly?

jmcastro2109 commented 2 years ago

Yes, it does!

Here is the code that is giving me the error

`

Trying Stochastic Programs

using JuMP, LinearAlgebra, Random, Distributions, BenchmarkTools, MosekTools, NLopt

using StochasticPrograms

@sampler SimpleSampler = begin N::MvNormal

SimpleSampler(μ, Σ) = new(MvNormal(μ, Σ))

@sample Scenario begin
    x = rand(sampler.N)
    return Scenario(z₁ = exp(x[1]), z₂ = exp(x[2]))
end

end

μ = [-0.5, -0.5] Σ = [1.0 0.0 0.0 1]

sampler = SimpleSampler(μ, Σ)

@stochastic_model simple_model begin @stage 1 begin @decision(simple_model, x₁ >= 0) @decision(simple_model, x₂ >= 0) @objective(simple_model, Max, -1x₁ - 1x₂) end @stage 2 begin @uncertain z₁ z₂ @recourse(simple_model, 0 <= ν₁ ) @recourse(simple_model, 0 <= ν₂ ) @recourse(simple_model, y[1:J]>= 0) @recourse(simple_model, t[1:J]>= 0) @objective(simple_model, Min, ν₁x₁+ ν₂x₂ + (((σ-1)/σ)^(σ-1)/(σ).β.^(σ))'y) for j=1:J @constraint(simple_model, τ[j,1]/z₁ + τ[j,1]ν₁ - t[j] >= 0) @constraint(simple_model, τ[j,2]/z₂ + τ[j,2]ν₂ - t[j] >= 0) end

    for j = 1 : J
        @constraint(simple_model, [y[j],t[j],1] in MOI.PowerCone(1/σ))
    end 
end

end

import Ipopt

sampled_sp = instantiate(simple_model, sampler, 5, optimizer = Mosek.Optimizer)

print(sampled_sp)

StochasticPrograms.optimize!(sampled_sp) `

martinbiel commented 2 years ago

Parameter values are missing.

It looks like the objective is not being bridged properly. Quadratic support is not thoroughly tested, and conic stuff are not tested at all, so it is not unlikely that this is the issue.