Help with formulation - Githubissues

Hi Martin, This is an amazing package for stochastic programming in Julia and I really like it. Thanks for the great contribution!

I have coded a CVaR optimisation problem using the extensive form (deterministic equivalent), but I would like to solve it using Lshaped method using your package. This is because I run into memory issues if number of scenarios or number of investments are increased.

I have given both my implementation using the the extensive form and my attempt at using your package. Would you be able to suggest corrections to the code in my attempt at using your package?

#Extensive form without using StochasticPrograms.jl
using JuMP
using Clp
# Generate random data for initial value of 5 investments
InitialValue= float(rand(80:140,1,5))
# Generate 10,000 scenarios for future value of investments and uncertain demand
FutureValue = float(rand(40:150,10000,5))
Demand = float(rand(40:50,10000))
nscen = size(FutureValue,1)
ninvestments = size(FutureValue,2)
# Create a vector of probabilities for scenarios 
prob = 1/nscen*ones(nscen)
# Calculate loss amount based on changes in value of investments
Loss = zeros(Float64,nscen,ninvestments)
Loss = -FutureValue.+InitialValue
# Assume 95% confidence level
α = 0.05
cvarRanmodel = Model(Clp.Optimizer)
@variable(cvarRanmodel,x[b=1:ninvestments]>=0)
@variable(cvarRanmodel,y[b=1:nscen]>=0)
@variable(cvarRanmodel,γ)
@constraint(cvarRanmodel,budget,sum(x[i] for i =1:ninvestments) == 1)
@constraint(cvarRanmodel,constraint2[j in 1:nscen], y[j]-sum(Loss[j,i]*x[i] for i in 1:ninvestments) + γ >= 0)
@constraint(cvarRanmodel,constraint3[j in 1:nscen], sum(FutureValue[j,i]*x[i] for i in 1:ninvestments)-Demand[j] >= 0)
@objective(cvarRanmodel,Min,γ+1/(α)*sum(y[j]*prob[j] for j =1:nscen))
optimize!(cvarRanmodel)
println("")
println("VaR ", JuMP.value.(γ))
println("CVaR ", objective_value(cvarRanmodel))
println("Investment allocation",JuMP.value.(x))

My attempt at formulating it using your package is given below:

using StochasticPrograms
using GLPK
# Generate random data for initial value of 5 investments
InitialValue= float(rand(80:140,1,5))
# Generate 10,000 scenarios for future value of investments and uncertain demand
FutureValue = float(rand(40:150,10000,5))
Demand = float(rand(40:50,10000))
nscen = size(FutureValue,1)
ninvestments = size(FutureValue,2)
# Create a vector of probabilities for scenarios 
prob = 1/nscen*ones(nscen)
# Calculate loss amount based on changes in value of investments
Loss = zeros(Float64,nscen,ninvestments)
Loss = -FutureValue.+InitialValue
# Assume 95% confidence level
α = 0.05

sp_model = @stochastic_model begin
    @stage 1 begin
        @decision(model, x[b=1:ninvestments]>=0)
        @variable(model,γ)
        @objective(model, Min, γ)
        @constraint(model, sum(x[i] for i =1:ninvestments) == 1)
    end
    @stage 2 begin
        @uncertain Loss FutureValue Demand
        @variable(model,y[b=1:nscen]>=0)
        @objective(model,Min,1/(α)*sum(y[j]*prob[j] for j =1:nscen)) #Not sure how to write in StochasticPrograms.jl format
        @constraint(model,constraint2[j in 1:nscen], y[j]-sum(Loss[j,i]*x[i] for i in 1:ninvestments) + γ >= 0)  #Not sure how to write in StochasticPrograms.jl format
        @constraint(model,constraint3[j in 1:nscen], sum(FutureValue[j,i]*x[i] for i in 1:ninvestments)-Demand[j] >= 0)  #Not sure how to write in StochasticPrograms.jl format
    end
end
ξ = Scenario(Loss,FutureValue, Demand, y,prob ) #This is the part I am not getting right
sp_lshaped = instantiate(model, [ξ], optimizer = LShaped.Optimizer)
set_optimizer_attribute(sp_lshaped, MasterOptimizer(), GLPK.Optimizer)
set_optimizer_attribute(sp_lshaped, SubproblemOptimizer(), GLPK.Optimizer)
optimize!(sp_lshaped)

martinbiel / StochasticPrograms.jl

Help with formulation #8