Skip to content

Commit e64dac0

Browse files
Henri-GerardHenri-Gerard
authored
Merge pull request #151 from JuliaOpt/add-risk
Add possibility to compute risk
2 parents 245986d + a7d43b4 commit e64dac0

File tree

11 files changed

+584
-13
lines changed

11 files changed

+584
-13
lines changed

‎examples/risk-example.jl‎

Lines changed: 162 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard
2+
# This Source Code Form is subject to the terms of the Mozilla Public
3+
# License, v. 2.0. If a copy of the MPL was not distributed with this
4+
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
5+
#############################################################################
6+
# Test impact of risk solving a stock problem :
7+
# Min F [\sum_{t=1}^TF c_t u_t]
8+
# s.t. s_{t+1} = s_t + u_t - xi_t, s_0 given
9+
# 0 <= s_t <= 1
10+
# u_min <= u_t <= u_max
11+
# u_t choosen knowing xi_1 .. xi_t
12+
#############################################################################
13+
14+
using StochDynamicProgramming, Clp
15+
println("library loaded")
16+
17+
run_Expectation = true # false if you don't want to test Expectation
18+
run_AVaR = true # false if you don't want to test AVaR
19+
run_WorstCase = true # false if you don't want to test WorstCase
20+
run_ConvexCombi = true # false if you don't want to test Convex Combination
21+
run_Polyhedral = true # false if you don't want to test Polyhedral Risk Measure
22+
23+
######## Optimization parameters ########
24+
# choose the LP solver used.
25+
SOLVER = ClpSolver() # require "using Clp"
26+
#const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0) # require "using CPLEX"
27+
28+
# convergence test
29+
MAX_ITER = 10 # number of iterations of SDDP
30+
31+
######## Stochastic Model Parameters ########
32+
N_STAGES = 6 # number of stages of the SP problem
33+
COSTS = [sin(3*t)-1 for t in 1:N_STAGES-1]
34+
#const COSTS = rand(N_STAGES) # randomly generating deterministic costs
35+
36+
CONTROL_MAX = 0.5 # bounds on the control
37+
CONTROL_MIN = 0
38+
39+
XI_MAX = 0.3 # bounds on the noise
40+
XI_MIN = 0
41+
N_XI = 10 # discretization of the noise
42+
43+
S0 = 0.5 # initial stock
44+
45+
# create law of noises
46+
proba = 1/N_XI*ones(N_XI) # uniform probabilities
47+
xi_support = collect(linspace(XI_MIN,XI_MAX,N_XI))
48+
xi_law = NoiseLaw(xi_support, proba)
49+
xi_laws = NoiseLaw[xi_law for t in 1:N_STAGES-1]
50+
51+
# Define dynamic of the stock:
52+
function dynamic(t, x, u, xi)
53+
return [x[1] + u[1] - xi[1]]
54+
end
55+
56+
# Define cost corresponding to each timestep:
57+
function cost_t(t, x, u, w)
58+
return COSTS[t] * u[1]
59+
end
60+
61+
######## Setting up the SPmodel
62+
s_bounds = [(0, 1)] # bounds on the state
63+
u_bounds = [(CONTROL_MIN, CONTROL_MAX)] # bounds on controls
64+
65+
println("Initializing functions to compare execution time")
66+
spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = Expectation())
67+
set_state_bounds(spmodel, s_bounds) # adding the bounds to the model
68+
# 10 forward pass, stop at MAX_ITER
69+
paramSDDP = SDDPparameters(SOLVER,
70+
passnumber=10,
71+
max_iterations=MAX_ITER)
72+
sddp = solve_SDDP(spmodel, paramSDDP, 0) # display information every 2 iterations
73+
lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
74+
75+
######### Solving the problem via SDDP with Expectation
76+
if run_Expectation
77+
tic()
78+
spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = Expectation())
79+
set_state_bounds(spmodel, s_bounds) # adding the bounds to the model
80+
println("Expectation's model set up")
81+
println("Starting resolution with Expectation")
82+
# 10 forward pass, stop at MAX_ITER
83+
paramSDDP = SDDPparameters(SOLVER,
84+
passnumber=10,
85+
max_iterations=MAX_ITER)
86+
sddp = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
87+
lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
88+
println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
89+
toc(); println();
90+
end
91+
92+
######### Solving the problem via SDDP with AVaR
93+
if run_AVaR
94+
tic()
95+
spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = AVaR((N_XI-1)/N_XI))
96+
set_state_bounds(spmodel, s_bounds) # adding the bounds to the model
97+
println("AVaR's model set up")
98+
println("Starting resolution with AVaR")
99+
# 10 forward pass, stop at MAX_ITER
100+
paramSDDP = SDDPparameters(SOLVER,
101+
passnumber=10,
102+
max_iterations=MAX_ITER)
103+
sddp = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
104+
lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
105+
println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
106+
toc(); println();
107+
end
108+
109+
######### Solving the problem via SDDP with Worst Case
110+
if run_WorstCase
111+
tic()
112+
spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = WorstCase())
113+
set_state_bounds(spmodel, s_bounds) # adding the bounds to the model
114+
println("Worst Case's model set up")
115+
println("Starting resolution with Worst Case")
116+
# 10 forward pass, stop at MAX_ITER
117+
paramSDDP = SDDPparameters(SOLVER,
118+
passnumber=10,
119+
max_iterations=MAX_ITER)
120+
sddp = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
121+
lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
122+
println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
123+
toc(); println();
124+
end
125+
126+
if run_ConvexCombi
127+
tic()
128+
spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = ConvexCombi((N_XI-1)/N_XI,0.5))
129+
set_state_bounds(spmodel, s_bounds) # adding the bounds to the model
130+
println("Convex Combination's model set up")
131+
println("Starting resolution with Convex Combination")
132+
# 10 forward pass, stop at MAX_ITER
133+
paramSDDP = SDDPparameters(SOLVER,
134+
passnumber=10,
135+
max_iterations=MAX_ITER)
136+
sddp = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
137+
lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
138+
println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
139+
toc(); println();
140+
end
141+
142+
if run_Polyhedral
143+
beta = (N_XI-1)/N_XI
144+
prob = 1/N_XI*ones(N_XI)
145+
polyset = repmat(1/beta*prob',N_XI)
146+
for i = 1:N_XI
147+
polyset[i,i] = (beta*N_XI-N_XI+1)/(N_XI*beta)
148+
end
149+
tic()
150+
spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = PolyhedralRisk(polyset))
151+
set_state_bounds(spmodel, s_bounds) # adding the bounds to the model
152+
println("Polyhedral's model set up")
153+
println("Starting resolution with Polyhedral")
154+
# 10 forward pass, stop at MAX_ITER
155+
paramSDDP = SDDPparameters(SOLVER,
156+
passnumber=10,
157+
max_iterations=MAX_ITER)
158+
sddp = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
159+
lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
160+
println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
161+
toc(); println();
162+
end

‎src/SDDPoptimize.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ function solve_SDDP(model::SPModel, param::SDDPparameters, verbosity=0::Int64, v
4040
stopcrit::AbstractStoppingCriterion=IterLimit(param.max_iterations),
4141
prunalgo::AbstractCutPruningAlgo=CutPruners.AvgCutPruningAlgo(-1),
4242
regularization=nothing)
43+
4344
sddp = SDDPInterface(model, param,
4445
stopcrit,
4546
prunalgo,

‎src/StochDynamicProgramming.jl‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,9 @@ export solve_SDDP,
2424
PolyhedralFunction, forward_simulations,
2525
StochDynProgModel, SDPparameters, solve_dp,
2626
sampling, get_control, get_bellman_value,
27-
benchmark_parameters, SDDPInterface
27+
benchmark_parameters, SDDPInterface,
28+
risk_proba,
29+
RiskMeasure, AVaR, Expectation, WorstCase, ConvexCombi, PolyhedralRisk
2830

2931
include("noises.jl")
3032
include("objects.jl")
@@ -41,4 +43,5 @@ include("sdp.jl")
4143
include("compare.jl")
4244
include("cutpruning.jl")
4345
include("simulation.jl")
46+
include("risk.jl")
4447
end

‎src/forwardBackwardIterations.jl‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -303,7 +303,6 @@ function backward_pass!(sddp::SDDPInterface,
303303
sddp.stats.solverexectime_bw = vcat(sddp.stats.solverexectime_bw, solvertime)
304304
end
305305

306-
307306
"""Compute cuts in Hazard-Decision (classical SDDP)."""
308307
function compute_cuts_hd!(model::SPModel, param::SDDPparameters,
309308
V::Vector{PolyhedralFunction},
@@ -345,6 +344,9 @@ function compute_cuts_hd!(model::SPModel, param::SDDPparameters,
345344
# Scale probability (useful when some problems where infeasible):
346345
proba /= sum(proba)
347346

347+
#Modify the probability vector to compute the value of the risk measure
348+
proba = risk_proba(proba,model.riskMeasure,costs)
349+
348350
# Compute expectation of subgradient λ:
349351
subgradient = vec(sum(proba' .* subgradient_array, 2))
350352
# ... expectation of cost:

‎src/objects.jl‎

Lines changed: 66 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,57 @@
77
#############################################################################
88

99

10+
@compat abstract type RiskMeasure end
11+
12+
# Define an object to
13+
type Expectation <: RiskMeasure
14+
function Expectation()
15+
return new()
16+
end
17+
end
18+
19+
type AVaR <: RiskMeasure
20+
# If the random variable is a cost and beta = 0.05,
21+
# it returns the average of the five worst costs solving the problem
22+
# minimize alpha + 1/beta * E[max(X - alpha; 0)]
23+
# beta = 1 ==> AVaR = Expectation
24+
# beta = 0 ==> AVaR = WorstCase
25+
beta::Float64
26+
function AVaR(beta)
27+
return new(beta)
28+
end
29+
end
30+
31+
type WorstCase <: RiskMeasure
32+
function WorstCase()
33+
return new()
34+
end
35+
end
36+
37+
type ConvexCombi <: RiskMeasure
38+
# Define a convex combination between Expectation and AVaR_{beta}
39+
# with form lambda*E + (1-lambda)*AVaR
40+
# lambda = 1 ==> Expectation
41+
# lambda = 0 ==> AVaR
42+
beta::Float64
43+
lambda::Float64
44+
function ConvexCombi(beta,lambda)
45+
return new(beta,lambda)
46+
end
47+
end
48+
49+
type PolyhedralRisk <: RiskMeasure
50+
# Define a convex polyhedral set P of probability distributions
51+
# by its extreme points p1, ..., pn
52+
# In the case of costs X, the problem solved is
53+
# maximize E_{pi}[X]
54+
# p1, ..., pn
55+
polyset::Array{Float64,2}
56+
function PolyhedralRisk(polyset)
57+
return new(polyset)
58+
end
59+
end
60+
1061
@compat abstract type SPModel end
1162

1263

@@ -57,17 +108,22 @@ type LinearSPModel <: SPModel
57108

58109
IS_SMIP::Bool
59110

111+
#Define the risk measure used at each stage
112+
riskMeasure::RiskMeasure
113+
60114
function LinearSPModel(n_stage, # number of stages
61115
u_bounds, # bounds of control
62-
x0, # initial state
63-
cost, # cost function
64-
dynamic, # dynamic
65-
aleas; # modelling of noises
66-
Vfinal=nothing, # final cost
67-
eqconstr=nothing, # equality constraints
68-
ineqconstr=nothing, # inequality constraints
69-
info=:HD, # information structure
70-
control_cat=nothing) # category of controls
116+
x0, # initial state
117+
cost, # cost function
118+
dynamic, # dynamic
119+
aleas; # modelling of noises
120+
Vfinal=nothing, # final cost
121+
eqconstr=nothing, # equality constraints
122+
ineqconstr=nothing, # inequality constraints
123+
info=:HD, # information structure
124+
control_cat=nothing, # category of controls
125+
riskMeasure = Expectation()
126+
)
71127

72128
# infer the problem's dimension
73129
dimStates = length(x0)
@@ -89,7 +145,7 @@ type LinearSPModel <: SPModel
89145
x_bounds = [(-Inf, Inf) for i=1:dimStates]
90146

91147
return new(n_stage, dimControls, dimStates, dimNoises, x_bounds, u_bounds,
92-
x0, cost, dynamic, aleas, Vf, isbu, eqconstr, ineqconstr, info, is_smip)
148+
x0, cost, dynamic, aleas, Vf, isbu, eqconstr, ineqconstr, info, is_smip, riskMeasure)
93149
end
94150
end
95151

0 commit comments

Comments
 (0)