[DEV] Take comments of pull request into account

Author: Henri-Gerard

examples/risk-example.jl

@@ -61,7 +61,7 @@ s_bounds = [(0, 1)] 			# bounds on the state
 u_bounds = [(CONTROL_MIN, CONTROL_MAX)] # bounds on controls
 
 println("Initialzing functions to compare execution time")
-spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, Expectation())
+spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = Expectation())
 set_state_bounds(spmodel, s_bounds) 	# adding the bounds to the model
 # 10 forward pass, stop at MAX_ITER
 paramSDDP = SDDPparameters(SOLVER,
@@ -73,7 +73,7 @@ lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellm
 ######### Solving the problem via SDDP with Expectation
 if run_expectation
     tic()
-    spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, Expectation())
+    spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = Expectation())
     set_state_bounds(spmodel, s_bounds) 	# adding the bounds to the model
     println("Expectation's model set up")
     println("Starting resolution with Expectation")
@@ -107,7 +107,7 @@ end
 ######### Solving the problem via SDDP with Worst Case
 if run_WorstCase
     tic()
-    spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, WorstCase())
+    spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = WorstCase())
     set_state_bounds(spmodel, s_bounds) 	# adding the bounds to the model
     println("Worst Case's model set up")
     println("Starting resolution with Worst Case")

examples/stock-example.jl

@@ -61,7 +61,7 @@ end
 ######## Setting up the SPmodel
 s_bounds = [(0, 1)] 			# bounds on the state
 u_bounds = [(CONTROL_MIN, CONTROL_MAX)] # bounds on controls
-spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, Expectation())
+spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws)
 set_state_bounds(spmodel, s_bounds) 	# adding the bounds to the model
 println("Model set up")
 

src/SDDPoptimize.jl

@@ -40,6 +40,7 @@ function solve_SDDP(model::SPModel, param::SDDPparameters, verbosity=0::Int64, v
                     stopcrit::AbstractStoppingCriterion=IterLimit(param.max_iterations),
                     prunalgo::AbstractCutPruningAlgo=CutPruners.AvgCutPruningAlgo(-1),
                     regularization=nothing)
+                    
     sddp = SDDPInterface(model, param,
                          stopcrit,
                          prunalgo,

src/StochDynamicProgramming.jl

@@ -25,8 +25,8 @@ export solve_SDDP,
         StochDynProgModel, SDPparameters, solve_dp,
         sampling, get_control, get_bellman_value,
         benchmark_parameters, SDDPInterface,
-        change_proba_risk,
-        RiskMeasure, CVaR, Expectation, WorstCase, ConvexCombi
+        argsup_proba_risk,
+        RiskMeasure, AVaR, Expectation, WorstCase, ConvexCombi
 
 include("noises.jl")
 include("objects.jl")

src/forwardBackwardIterations.jl

@@ -345,7 +345,7 @@ function compute_cuts_hd!(model::SPModel, param::SDDPparameters,
         proba /= sum(proba)
 
         #Modify the probability vector to compute the value of the risk measure
-        proba = change_proba_risk(proba,model.riskMeasure,sortperm(costs,rev = true))
+        proba = argsup_proba_risk(proba,model.riskMeasure,costs)
 
         # Compute expectation of subgradient λ:
         subgradient = vec(sum(proba' .* subgradient_array, 2))

src/objects.jl

@@ -16,10 +16,14 @@ type Expectation <: RiskMeasure
     end
 end
 
-type CVaR <: RiskMeasure
-    # level of risk aversity (0 = Expectation, 1 = Worst Case):
+type AVaR <: RiskMeasure
+    # If the random variable is a cost and beta = 0.05,
+    # it returns the average of the five worst costs solving the problem
+    # minimize alpha + 1/beta * E[max(X - alpha; 0)]
+    # beta = 1 ==> AVaR = Expectation
+    # beta = 0 ==> AVaR = WorstCase
     beta::Float64
-    function CVaR(beta)
+    function AVaR(beta)
         return new(beta)
     end
 end
@@ -31,15 +35,28 @@ type WorstCase <: RiskMeasure
 end
 
 type ConvexCombi <: RiskMeasure
-    # level of risk aversity (0 = Expectation, 1 = Worst Case):
+    # Define a convex combination between Expectation and AVaR_{beta}
+    # with form lambda*E + (1-lambda)*AVaR
+    # lambda = 1 ==> Expectation
+    # lambda = 0 ==> AVaR
     beta::Float64
-    #Convex coefficient
     lambda::Float64
     function ConvexCombi(beta,lambda)
         return new(beta,lambda)
     end
 end
 
+type PolyhedralRisk <: RiskMeasure
+    # Define a convex polyhedral set P of probability distributions
+    # by its extreme points p1, ..., pn
+    # In the case of costs X, the problem solved is
+    #   maximize     E_{pi}[X]
+    # p1, ..., pn
+    polyset::Array{Float64,2}
+    function PolyhedralRisk(polyset)
+        return new(polyset)
+    end
+end
 
 @compat abstract type SPModel end
 
@@ -100,12 +117,12 @@ type LinearSPModel <: SPModel
                            cost,                # cost function
                            dynamic,             # dynamic
                            aleas,               # modelling of noises
-                           riskMeasure;
                            Vfinal=nothing,      # final cost
                            eqconstr=nothing,    # equality constraints
                            ineqconstr=nothing,  # inequality constraints
                            info=:HD,            # information structure
                            control_cat=nothing, # category of controls
+                           riskMeasure = Expectation();
                            )
 
         # infer the problem's dimension

src/risk.jl

@@ -7,38 +7,44 @@
 #############################################################################
 
 """
-Update probability distribution to consider only worst case scenarios
+Risk measures consider here can be written (in case of a minimization problem)
+   sup     E_{p}[x]
+  p in P
+where P is a convex set of probability distribution.
+The function returns the argsup allowing to compute the
+risk measure as a expectation
 
 $(SIGNATURES)
 
 # Description
-Keep only the tail of the probability distribution given an order
-and renormalized it
+Return the probability among the set of possible probabilities
+allowing to compute the risk
 
 # Arguments
-* `prob::probabilitydistribution`:
-    SDDP interface object
-* `beta::Float64`:
-    real number
-* `perm::Array{float,1}`:
-    array of permutations
+* `prob::Array{float,1}`:
+    probability distribution
+* `riskMeasure::RiskMeasure`:
+    risk Measure
+* `costs`:
+    array of costs
 
 # Returns
-* `aux::Array{float,1}`:
+* `proba::Array{float,1}`:
     a new probability distribution taking risk into account.
 
 """
-function change_proba_risk(prob, riskMeasure::RiskMeasure,perm)
-    erro("'change_proba_risk' not defined for $(typedof(s))")
+function argsup_proba_risk(prob, riskMeasure::RiskMeasure,costs)
+    error("'argsup_proba_risk' not defined for $(typedof(s))")
 end
 
 """
 $(TYPEDEF)
 
-Update the probability distribution to calculate a Conditional Value at Risk of level beta.
+Return the probability distribution to compute a Average Value at Risk of level beta.
 """
-function change_proba_risk(prob,riskMeasure::CVaR,perm)
-    beta = 1-riskMeasure.beta
+function argsup_proba_risk(prob,riskMeasure::AVaR,costs)
+    perm = sortperm(costs,rev = true)
+    beta = riskMeasure.beta
     proba = zeros(length(prob))
     for i in 1:length(perm)
         proba[i] = prob[perm[i]]
@@ -61,31 +67,32 @@ end
 """
 $(TYPEDEF)
 
-Leave the probability distribution unchanged to calculate the expectation.
+Leave the probability distribution unchanged to compute the expectation.
 """
-function change_proba_risk(prob,riskMeasure::Expectation,perm)
+function argsup_proba_risk(prob,riskMeasure::Expectation,costs)
     return prob
 end
 
 """
 $(TYPEDEF)
-
+perm
 Return a dirac on the worst cost as a probability distribution.
 """
-function change_proba_risk(prob,riskMeasure::WorstCase,perm)
+function argsup_proba_risk(prob,riskMeasure::WorstCase,costs)
     proba = zeros(length(prob))
-    proba[perm[1]] = 1
+    proba[indmax(prob)] = 1
     return proba
 end
 
 """
 $(TYPEDEF)
 
-Update the probability distribution to calculate a convex combination
-between expactation and a Condiational Value at Risk of level beta.
+Return the probability distribution to compute a convex combination
+between expactation and an Average Value at Risk of level beta.
 """
-function change_proba_risk(prob,riskMeasure::ConvexCombi,perm)
-    beta = 1-riskMeasure.beta
+function argsup_proba_risk(prob,riskMeasure::ConvexCombi,costs)
+    perm = sortperm(costs,rev = true)
+    beta = riskMeasure.beta
     lambda = riskMeasure.lambda
     proba = zeros(length(prob))
     for i in 1:length(perm)
@@ -105,3 +112,15 @@ function change_proba_risk(prob,riskMeasure::ConvexCombi,perm)
     end
     return lambda*prob + (1-lambda)*aux
 end
+
+"""
+$(TYPEDEF)
+
+Return the worst extreme probability distribution
+defining the convex set P
+"""
+function argsup_proba_risk(prob,riskMeasure::PolyhedralRisk,costs)
+    P = riskMeasure.polyset
+    valuesup = P*costs
+    return P[indmax(valuesup),:]
+end

src/stopcrit.jl

@@ -137,11 +137,11 @@ Stops if the lower bound is stabilized
 total time of execution is greater than the time limit specified.
 For instance, `TimeLimit(100)` stops after 100s.
 """
-type Stabilization <: AbstractStoppingCriterion
+type LBStabilization <: AbstractStoppingCriterion
     epsilon::Float64
     n_back::Int
 end
 
-function stop(s::Stabilization, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats)
+function stop(s::LBStabilization, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats)
     totalstats.niterations > s.n_back && ((stats.lowerbound[end] - stats.lowerbound[end-s.n_back]) < epsilon)
 end

test/changeprob.jl

@@ -3,54 +3,64 @@ using Base.Test
 
 EPSILON = 0.0001
 
+# Check that there is no problem of definition
 @testset "Risk measure definition" begin
     @test isa(Expectation(), RiskMeasure)
     @test isa(WorstCase(), RiskMeasure)
-    for i in 0:10
-        @test isa(CVaR(i/10), RiskMeasure)
-        for j in 0:10
-          @test isa(ConvexCombi(i/10, j/10), RiskMeasure)
-        end
-    end
+    @test isa(AVaR(0), RiskMeasure)
+    @test isa(AVaR(0.5), RiskMeasure)
+    @test isa(AVaR(1), RiskMeasure)
+    @test isa(ConvexCombi(0,0), RiskMeasure)
+    @test isa(ConvexCombi(0,1), RiskMeasure)
+    @test isa(ConvexCombi(1,0), RiskMeasure)
+    @test isa(ConvexCombi(1,1), RiskMeasure)
+    @test isa(ConvexCombi(0.5,0.5), RiskMeasure)
 end
 
-
+# Test limit cases
+# An AVaR with beta = 0 is a WorstCase
+# An AVaR with beta = 1 is an Expectation
+# A ConvexCombi with lambda = 0 is an AVaR
+# A ConvexCombi with lambda = 1 is an Expectation
 @testset "Equality formulations" begin
-    @testset "Equality WorstCase CVaR(1)" begin
+    @testset "Equality WorstCase AVaR(0)" begin
         n = 100
         prob = 1/n*ones(n)
-        @test sum(abs.(change_proba_risk(prob,CVaR(1),1:n)-change_proba_risk(prob,WorstCase(),1:n)) .<= EPSILON*ones(n)) == n
+        @test sum(abs.(argsup_proba_risk(prob,AVaR(0),1:n)-argsup_proba_risk(prob,WorstCase(),1:n)) .<= EPSILON*ones(n)) == n
     end
 
-    @testset "Equality Expectation CVaR(0)" begin
+    @testset "Equality Expectation AVaR(1)" begin
         n = 100
         prob = 1/n*ones(n)
-        @test sum(abs.(change_proba_risk(prob,CVaR(0),1:n)-change_proba_risk(prob,Expectation(),1:n)) .<= EPSILON*ones(n)) == n
+        @test sum(abs.(argsup_proba_risk(prob,AVaR(0),1:n)-argsup_proba_risk(prob,Expectation(),1:n)) .<= EPSILON*ones(n)) == n
     end
 
     @testset "Equality Expectation ConvexCombi(beta,1)" begin
         n = 100
         prob = 1/n*ones(n)
-        @test sum(abs.(change_proba_risk(prob,Expectation(),1:n)-change_proba_risk(prob,ConvexCombi(rand(),1),1:n)) .<= EPSILON*ones(n)) == n
+        @test sum(abs.(argsup_proba_risk(prob,Expectation(),1:n)-argsup_proba_risk(prob,ConvexCombi(rand(),1),1:n)) .<= EPSILON*ones(n)) == n
     end
 
-    @testset "Equality CVaR(beta) ConvexCombi(beta,0)" begin
+    @testset "Equality AVaR(beta) ConvexCombi(beta,0)" begin
         n = 100
         beta = rand()
         prob = 1/n*ones(n)
-        @test sum(abs.(change_proba_risk(prob,CVaR(beta),1:n)-change_proba_risk(prob,ConvexCombi(beta,0),1:n)) .<= EPSILON*ones(n)) == n
+        @test sum(abs.(argsup_proba_risk(prob,AVaR(beta),1:n)-argsup_proba_risk(prob,ConvexCombi(beta,0),1:n)) .<= EPSILON*ones(n)) == n
     end
 
-    @testset "Right sense of CVaR" begin
-        betamin = rand()/2
-        betamax = rand()/2+0.5
+    # Check that in the case of minimization, AVaR find the worst costs
+    @testset "Right sense of AVaR" begin
+        betamax = rand()/2
+        betamin = rand()/2+0.5
         n = 100
         prob = 1/n*ones(n)
-        @test change_proba_risk(prob, CVaR(betamax), 1:n)'*(n:-1:1) - change_proba_risk(prob, CVaR(betamin), 1:n)'*(n:-1:1) >= 0
+        @test argsup_proba_risk(prob, AVaR(betamax), 1:n)'*(n:-1:1) - argsup_proba_risk(prob, AVaR(betamin), 1:n)'*(n:-1:1) >= 0
     end
 end
 
-@testset "Equality CVaR linear program" begin
+# AVaR can be computed as a linear program
+# We check equality between the formulation of Rockafellar and Urysev
+@testset "Equality AVaR linear program" begin
     n = 100
     X = rand(-100:100,100)
     beta = rand()
@@ -63,11 +73,11 @@ end
 
     @constraint(m, theta[1:n] .>= X[1:n] - alpha)
 
-    @objective(m, Min, alpha + 1/(1-beta)*sum(prob[i]*theta[i] for i in 1:n))
+    @objective(m, Min, alpha + 1/(beta)*sum(prob[i]*theta[i] for i in 1:n))
 
     status = solve(m)
 
-    probaCVaR = change_proba_risk(prob, CVaR(beta), sortperm(X, rev = true))
+    probaAVaR = argsup_proba_risk(prob, AVaR(beta), sortperm(X, rev = true))
 
-    @test abs(probaCVaR'*X - getobjectivevalue(m)) <= EPSILON
+    @test abs(probaAVaR'*X - getobjectivevalue(m)) <= EPSILON
 end

test/extensive_formulation.jl

@@ -39,7 +39,7 @@ using StochDynamicProgramming, Base.Test, Clp
     u_bounds = [(0., 7.), (0., Inf)]
 
     model_ef = StochDynamicProgramming.LinearSPModel(n_stages, u_bounds,
-                                                        x0, cost, dynamic, laws, Expectation())
+                                                        x0, cost, dynamic, laws)
     x_bounds_ef = [(0, 100)]
     set_state_bounds(model_ef, x_bounds_ef)