[DEV] Add polyhedral risk measure

Author: Henri-Gerard

examples/risk-example.jl

@@ -14,9 +14,11 @@
 using StochDynamicProgramming, Clp
 println("library loaded")
 
-run_expectation = true # false if you don't want to test expectation
-run_AVaR        = true # false if you don't want to test AVaR
-run_WorstCase   = true # false if you don't want to test WorstCase
+run_Expectation   = true # false if you don't want to test Expectation
+run_AVaR          = true # false if you don't want to test AVaR
+run_WorstCase     = true # false if you don't want to test WorstCase
+run_ConvexCombi   = true # false if you don't want to test Convex Combination
+run_Polyhedral    = true # false if you don't want to test Polyhedral Risk Measure
 
 ######## Optimization parameters  ########
 # choose the LP solver used.
@@ -71,7 +73,7 @@ sddp = solve_SDDP(spmodel, paramSDDP, 0) # display information every 2 iteration
 lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
 
 ######### Solving the problem via SDDP with Expectation
-if run_expectation
+if run_Expectation
     tic()
     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
@@ -90,7 +92,7 @@ end
 ######### Solving the problem via SDDP with AVaR
 if run_AVaR
     tic()
-    spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = AVaR(0.5))
+    spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = AVaR((N_XI-1)/N_XI))
     set_state_bounds(spmodel, s_bounds) 	# adding the bounds to the model
     println("AVaR's model set up")
     println("Starting resolution with AVaR")
@@ -120,3 +122,41 @@ if run_WorstCase
     println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
     toc(); println();
 end
+
+if run_ConvexCombi
+    tic()
+    spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = ConvexCombi((N_XI-1)/N_XI,0.5))
+    set_state_bounds(spmodel, s_bounds) 	# adding the bounds to the model
+    println("Convex Combination's model set up")
+    println("Starting resolution with Convex Combination")
+    # 10 forward pass, stop at MAX_ITER
+    paramSDDP = SDDPparameters(SOLVER,
+                               passnumber=10,
+                               max_iterations=MAX_ITER)
+    sddp = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
+    lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
+    println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
+    toc(); println();
+end
+
+if run_Polyhedral
+    beta = (N_XI-1)/N_XI
+    prob = 1/N_XI*ones(N_XI)
+    polyset = repmat(1/beta*prob',N_XI)
+    for i = 1:N_XI
+        polyset[i,i] = (beta*N_XI-N_XI+1)/(N_XI*beta)
+    end
+    tic()
+    spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws, riskMeasure = PolyhedralRisk(polyset))
+    set_state_bounds(spmodel, s_bounds) 	# adding the bounds to the model
+    println("Polyhedral's model set up")
+    println("Starting resolution with Polyhedral")
+    # 10 forward pass, stop at MAX_ITER
+    paramSDDP = SDDPparameters(SOLVER,
+                               passnumber=10,
+                               max_iterations=MAX_ITER)
+    sddp = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
+    lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions)
+    println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
+    toc(); println();
+end

src/StochDynamicProgramming.jl

@@ -26,7 +26,7 @@ export solve_SDDP,
         sampling, get_control, get_bellman_value,
         benchmark_parameters, SDDPInterface,
         argsup_proba_risk,
-        RiskMeasure, AVaR, Expectation, WorstCase, ConvexCombi
+        RiskMeasure, AVaR, Expectation, WorstCase, ConvexCombi, PolyhedralRisk
 
 include("noises.jl")
 include("objects.jl")

test/changeprob.jl

@@ -15,6 +15,7 @@ EPSILON = 0.0001
     @test isa(ConvexCombi(1,0), RiskMeasure)
     @test isa(ConvexCombi(1,1), RiskMeasure)
     @test isa(ConvexCombi(0.5,0.5), RiskMeasure)
+    @test isa(PolyhedralRisk(1/100*ones(1,100)), RiskMeasure)
 end
 
 # Test limit cases
@@ -81,3 +82,23 @@ end
 
     @test abs(probaAVaR'*X - getobjectivevalue(m)) <= EPSILON
 end
+
+# The convex set of probabilty distributions of AVaR_{beta} is defined by
+# {Q | dQ/dP <= 1/beta}
+# We check equality between polyhedral formulation and AVaR formulation
+@testset "Equality AVaR Polyhedral" begin
+    n = 10
+    X = rand(-10:10,10)
+    beta = (n-1)/n+rand()*(1-(n-1)/n)
+    prob = 1/n*ones(n)
+
+    polyset = repmat(1/beta*prob',n)
+    for i = 1:n
+        polyset[i,i] = (beta*n-n+1)/(n*beta)
+    end
+
+    probaAVaR = argsup_proba_risk(prob, AVaR(beta), X)
+    probaPolyhedral = argsup_proba_risk(prob, PolyhedralRisk(polyset), X)
+
+    @test sum(abs.(probaAVaR-probaPolyhedral) .<= EPSILON*ones(n)) == n
+end

test/sddprisk.jl

@@ -51,26 +51,26 @@ using StochDynamicProgramming, JuMP, Clp, CutPruners
 
     V = nothing
     model = StochDynamicProgramming.LinearSPModel(n_stages, u_bounds,
-                                                  x0, cost, dynamic, laws)
-    model.beta = 0.0
-    println(model.beta)
+                                                  x0, cost, dynamic, laws, riskMeasure = WorstCase())
+
+
     set_state_bounds(model, x_bounds)
     # Test error if bounds are not well specified:
-    
+
 
     # Generate scenarios for forward simulations:
     noise_scenarios = simulate_scenarios(model.noises,param.forwardPassNumber)
 
     sddp_costs = 0
 
-    
+
         # Compute bellman functions with SDDP:
         sddp = solve_SDDP(model, param, 0)
-        
+
 
         V = sddp.bellmanfunctions
         # Test if the first subgradient has the same dimension as state:
-        
+
 
         # Test upper bounds estimation with Monte-Carlo:
         n_simulations = 100
@@ -79,29 +79,21 @@ using StochDynamicProgramming, JuMP, Clp, CutPruners
                                                            V,
                                                            sddp.solverinterface,
                                                            n_simulations)[1]
-        
+
 
         pbs = sddp.solverinterface
         sddp_costs, stocks = forward_simulations(model, param, pbs,
                                                  noise_scenarios)
         # Test error if scenarios are not given in the right shape:
-       
+
 
         # Test computation of optimal control:
         aleas = noise_scenarios[1, 1, :]
         opt = StochDynamicProgramming.get_control(model, param,
                                                   sddp.solverinterface,
                                                   1, model.initialState, aleas)
-        
+
 
         # Test display:
         param.compute_ub = 0
         sddp = solve_SDDP(model, param, V, 2)
-   
-
-
-
-
-
-
-