Merge pull request #107 from JuliaOpt/ef-fix

Ef fix

Author: leclere

examples/stock-example.jl

@@ -24,10 +24,10 @@ const SOLVER = ClpSolver()
 # const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0) # require "using CPLEX"
 
 # convergence test
-const MAX_ITER = 100 # number of iterations of SDDP
+const MAX_ITER = 10 # number of iterations of SDDP
 
 ######## Stochastic Model  Parameters  ########
-const N_STAGES = 5              # number of stages of the SP problem
+const N_STAGES = 6              # number of stages of the SP problem
 const COSTS = rand(N_STAGES)    # generating deterministic costs
 
 const CONTROL_MAX = 0.5         # bounds on the control
@@ -65,8 +65,8 @@ println("Model set up")
 ######### Solving the problem via SDDP
 if run_sddp
     println("Starting resolution by SDDP")
-    paramSDDP = SDDPparameters(SOLVER, 2, 0, MAX_ITER) # 2 forward pass, stop at MAX_ITER
-    V, pbs = solve_SDDP(spmodel, paramSDDP, 10) # display information every 10 iterations
+    paramSDDP = SDDPparameters(SOLVER, 10, 0, MAX_ITER) # 10 forward pass, stop at MAX_ITER
+    V, pbs = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations
     lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, V)
     println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))
 end
@@ -85,11 +85,10 @@ end
 
 if run_ef
     println("Starting resolution by Extensive Formulation")
-    println( extensive_formulation(spmodel, paramSDDP))
     value_ef = extensive_formulation(spmodel, paramSDDP)[1]
     println("Value obtained by Extensive Formulation: "*string(round(value_ef,4)))
-    println(round(value_sdp/value_ef,4))
-    println(round(lb_sddp/value_ef,4))
+    println("Relative error of SDP value: "*string(100*round(value_sdp/value_ef-1,4))*"%")
+    println("Relative error of SDDP lower bound: "*string(100*round(lb_sddp/value_ef-1,4))*"%")
 end
 
 ######### Comparing the solutions on simulated scenarios.
@@ -101,6 +100,6 @@ if run_sdp
     costsdp, states, stocks =sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs)
 end
 if run_sddp && run_sdp
-    println("Relative difference between sddp and sdp: "*string(2*round(mean(costsddp-costsdp)/mean(costsddp+costsdp),4)))
+    println("Simulated relative difference between sddp and sdp: "
+            *string(round(200*mean(costsddp-costsdp)/mean(costsddp+costsdp),3))*"%")
 end
-

src/SDDPoptimize.jl

@@ -39,6 +39,7 @@ fulfilled.
 function solve_SDDP(model::SPModel, param::SDDPparameters, display=0::Int64)
     # initialize value functions:
     V, problems = initialize_value_functions(model, param)
+    println("initial value function intialized")
     # Run SDDP upon example:
     sddp_stats = run_SDDP!(model, param, V, problems, display)
     return V, problems, sddp_stats
@@ -128,10 +129,16 @@ function run_SDDP!(model::SPModel,
         push!(stats.upper_bounds, upb)
 
         if (display > 0) && (iteration_count%display==0)
-            println("Pass number ", iteration_count,
+            if upb < Inf
+                println("Pass number ", iteration_count,
                     "\tUpper-bound: ", upb,
                     "\tLower-bound: ", round(stats.lower_bounds[end], 4),
                     "\tTime: ", round(stats.exectime[end], 2),"s")
+             else 
+                println("Pass number ", iteration_count,
+                    "\tLower-bound: ", round(stats.lower_bounds[end], 4),
+                    "\tTime: ", round(stats.exectime[end], 2),"s")
+             end
         end
 
     end
@@ -313,7 +320,7 @@ function initialize_value_functions(model::SPModel,
                                     param::SDDPparameters)
 
     solverProblems = build_models(model, param)
-
+    println("model builded")
     V = PolyhedralFunction[
                 PolyhedralFunction([], Array{Float64}(0, model.dimStates), 0) for i in 1:model.stageNumber]
 

src/extensiveFormulation.jl

@@ -53,12 +53,16 @@ function extensive_formulation(model, param)
             end
         end
     end
-
-    #Instantiate the problem creating dynamic constraint at each node
-    for t = 1 : (T)
+    #Add state constraints
+    for t = 1 :(T+1)
         for n = 1 : N[t]
             @constraint(mod,[x[t,DIM_STATE*(n-1)+k] for k = 1:DIM_STATE] .>= [model.xlim[k][1] for k = 1:DIM_STATE])
             @constraint(mod,[x[t,DIM_STATE*(n-1)+k] for k = 1:DIM_STATE] .<= [model.xlim[k][2] for k = 1:DIM_STATE])
+        end
+    end
+    #Instantiate the problem creating dynamic constraint at each node
+    for t = 1 : (T)
+        for n = 1 : N[t]
             for xi = 1 : laws[t].supportSize
                 m = (n-1)*laws[t].supportSize+xi
 
@@ -87,12 +91,18 @@ function extensive_formulation(model, param)
     @constraint(mod, [x[1,k] for k = 1:DIM_STATE] .== X_init)
 
     #Define the objective of the function
-    @objective(mod, Min, sum{ sum{proba[t][laws[t].supportSize*(n-1)+k]*c[t,laws[t].supportSize*(n-1)+k],k = 1:laws[t].supportSize} , t = 1:T, n=1:div(N[t+1],laws[t].supportSize)})
+    @objective(mod, Min, 
+    sum{ 
+        sum{    proba[t][laws[t].supportSize*(n-1)+k]*c[t,laws[t].supportSize*(n-1)+k],
+            k = 1:laws[t].supportSize}, 
+        t = 1:T, n=1:div(N[t+1],laws[t].supportSize)}
+    )
 
     status = solve(mod)
     solved = (status == :Optimal)
 
     if solved
+        println("EF value: "*string(getobjectivevalue(mod)))
         firstControl = collect(values(getvalue(u)))[1:DIM_CONTROL*laws[1].supportSize]
         return getobjectivevalue(mod), firstControl, status
     else