Merge pull request #117 from JuliaOpt/quad-regularization-repack

Quad regularization repack
+ add quadratic regularization (following T.Asamov work)
+ refactor run_sddp to better identify each step of the algorithm

Author: leclere

src/SDDPoptimize.jl

@@ -32,14 +32,11 @@ fulfilled.
 * `problems::Array{JuMP.Model}`:
     the collection of linear problems used to approximate
     each value function
-* `count_callsolver::Int64`:
-    number of times the solver has been called
+* `sddp_stats::SDDPStat`:
 
 """
 function solve_SDDP(model::SPModel, param::SDDPparameters, verbose=0::Int64)
-    if model.IS_SMIP && isa(param.MIPSOLVER, Void)
-        error("MIP solver is not defined. Please set `param.MIPSOLVER`")
-    end
+    check_SDDPparameters(model,param,verbose)
     # initialize value functions:
     V, problems = initialize_value_functions(model, param)
     (verbose > 0) && println("Initial value function loaded into memory.")
@@ -48,164 +45,145 @@ function solve_SDDP(model::SPModel, param::SDDPparameters, verbose=0::Int64)
     return V, problems, sddp_stats
 end
 
+"""
+Solve SDDP algorithm with hotstart and return estimation of bellman functions.
+
+# Description
+Alternate forward and backward phase till the stopping criterion is
+fulfilled.
 
+# Arguments
+* `model::SPmodel`:
+    the stochastic problem we want to optimize
+* `param::SDDPparameters`:
+    the parameters of the SDDP algorithm
+* `V::Vector{PolyhedralFunction}`:
+    current lower approximation of Bellman functions
+* `verbose::Int64`:
+    Default is `0`
+    If non null, display progression in terminal every
+    `n` iterations, where `n` is the number specified by display.
+
+# Returns
+* `V::Array{PolyhedralFunction}`:
+    the collection of approximation of the bellman functions
+* `problems::Array{JuMP.Model}`:
+    the collection of linear problems used to approximate
+    each value function
+* `sddp_stats::SDDPStat`:
+"""
 function solve_SDDP(model::SPModel, param::SDDPparameters, V::Vector{PolyhedralFunction}, verbose=0::Int64)
-    if model.IS_SMIP && isa(param.MIPSOLVER, Void)
-        error("MIP solver is not defined. Please set `param.MIPSOLVER`")
-    end
+    check_SDDPparameters(model,param,verbose)
     # First step: process value functions if hotstart is called
     problems = hotstart_SDDP(model, param, V)
     sddp_stats = run_SDDP!(model, param, V, problems, verbose)
     return V, problems, sddp_stats
 end
 
 
-"""Run SDDP iterations."""
+"""Run SDDP iterations.
+
+# Arguments
+* `model::SPmodel`:
+    the stochastic problem we want to optimize
+* `param::SDDPparameters`:
+    the parameters of the SDDP algorithm
+* `V::Vector{PolyhedralFunction}`:
+    Polyhedral lower approximation of Bellman functions
+* `problems::Vector{JuMP.Model}`:
+* `verbose::Int64`:
+    Default is `0`
+    If non null, display progression in terminal every
+    `n` iterations, where `n` is the number specified by display.
+
+# Returns
+* `stats:SDDPStats`:
+    contains statistics of the current algorithm
+"""
 function run_SDDP!(model::SPModel,
                     param::SDDPparameters,
                     V::Vector{PolyhedralFunction},
                     problems::Vector{JuMP.Model},
                     verbose=0::Int64)
 
     #Initialization of the counter
-    stats = SDDPStat(0, [], [], [], 0)
+    stats = SDDPStat()
 
     (verbose > 0) && println("Initialize cuts")
 
     # If computation of upper-bound is needed, a set of scenarios is built
     # to keep always the same realization for upper bound estimation:
-    if param.compute_upper_bound > 0
-        upperbound_scenarios = simulate_scenarios(model.noises, param.monteCarloSize)
-    end
+    #if param.compute_ub > 0 #TODO
+    upperbound_scenarios = simulate_scenarios(model.noises, param.in_iter_mc)
 
     upb = Inf
     costs = nothing
     stopping_test::Bool = false
-    iteration_count::Int64 = 0
+
 
     # Launch execution of forward and backward passes:
-    while (iteration_count < param.maxItNumber) & (~stopping_test)
+    while (~stopping_test)
         # Time execution of current pass:
         tic()
 
-        # Draw a set of scenarios according to the probability
-        # law specified in model.noises:
-        noise_scenarios = simulate_scenarios(model.noises, param.forwardPassNumber)
-
         ####################
-        # Forward pass
-        _, stockTrajectories,_,callsolver_forward = forward_simulations(model,
-                            param,
-                            problems,
-                            noise_scenarios)
+        # Forward pass : compute stockTrajectories
+        costs, stockTrajectories, callsolver_forward = forward_pass!(model, param, V, problems)
 
         ####################
-        # Backward pass
-        callsolver_backward = backward_pass!(model,
-                      param,
-                      V,
-                      problems,
-                      stockTrajectories,
-                      model.noises)
-
-        # Update the number of solver call
-        stats.ncallsolver += callsolver_forward + callsolver_backward
-        iteration_count += 1
-        stats.niterations += 1
+        # Backward pass : update polyhedral approximation of Bellman functions
+        callsolver_backward = backward_pass!(model, param, V, problems, stockTrajectories, model.noises)
 
         ####################
-        # If specified, prune cuts:
-        if (param.compute_cuts_pruning > 0) && (iteration_count%param.compute_cuts_pruning==0)
-            (verbose > 0) && println("Prune cuts ...")
-            remove_redundant_cuts!(V)
-            prune_cuts!(model, param, V)
-            problems = hotstart_SDDP(model, param, V)
-        end
-
-        # Update estimation of lower bound:
-        lwb = get_bellman_value(model, param, 1, V[1], model.initialState)
-        push!(stats.lower_bounds, lwb)
+        # cut pruning
+        prune_cuts!(model,param,V,stats.niterations,verbose)
 
         ####################
-        # If specified, compute upper-bound:
-        if (param.compute_upper_bound > 0) && (iteration_count%param.compute_upper_bound==0)
-            (verbose > 0) && println("Compute upper-bound with ",
-                                      param.monteCarloSize, " scenarios...")
-            # estimate upper-bound with Monte-Carlo estimation:
-            upb, costs = estimate_upper_bound(model, param, upperbound_scenarios, problems)
-            if param.gap > 0.
-                stopping_test = test_stopping_criterion(lwb, upb, param.gap)
-            end
-        end
+        # In iteration upper bound estimation
+        upb = in_iteration_upb_estimation(model, param, stats.niterations, verbose,
+                                          upperbound_scenarios, upb, problems)
 
-        push!(stats.exectime, toq())
-        push!(stats.upper_bounds, upb)
+        ####################
+        # Update stats
+        lwb = get_bellman_value(model, param, 1, V[1], model.initialState)
+        updateSDDPStat!(stats, callsolver_forward+callsolver_backward, lwb, upb, toq())
 
-        if (verbose > 0) && (iteration_count%verbose==0)
-            print("Pass number ", iteration_count)
-            (upb < Inf) && print("\tUpper-bound: ", upb)
-            println("\tLower-bound: ", round(stats.lower_bounds[end], 4),
-                "\tTime: ", round(stats.exectime[end], 2),"s")
-        end
+        print_current_stats(stats,verbose)
 
+        ####################
+        # Stopping test
+        stopping_test = test_stopping_criterion(param,stats)
+        stats.niterations += 1
     end
 
     ##########
     # Estimate final upper bound with param.monteCarloSize simulations:
-    if (verbose>0) && (param.compute_upper_bound >= 0)
-        V0 = get_bellman_value(model, param, 1, V[1], model.initialState)
+    display_final_solution(model, param,V,problems,stats,verbose)
+    return stats
+end
+
 
-        if param.compute_upper_bound == 0
+"""Display final results once SDDP iterations are finished."""
+function display_final_solution(model::SPModel, param::SDDPparameters, V, problems, stats::SDDPStat, verbose::Int64)
+    if (verbose>0) && (param.compute_ub >= 0)
+        lwb = get_bellman_value(model, param, 1, V[1], model.initialState)
+
+        if param.compute_ub == 0
             println("Estimate upper-bound with Monte-Carlo ...")
-            upb, costs = estimate_upper_bound(model, param, V, problems)
+            upb, costs = estimate_upper_bound(model, param, V, problems, param.monteCarloSize)
+        else
+            upb = stats.upperbounds[end]
         end
 
         println("Estimation of upper-bound: ", round(upb,4),
-                "\tExact lower bound: ", round(V0,4),
-                "\t Gap <  ", round(100*(upb-V0)/V0, 2) , "\%  with prob. > 97.5 \%")
+                "\tExact lower bound: ", round(lwb,4),
+                "\t Gap <  ", round(100*(upb-lwb)/lwb, 2) , "\%  with prob. > 97.5 \%")
         println("Estimation of cost of the solution (fiability 95\%):",
                  round(mean(costs),4), " +/- ", round(1.96*std(costs)/sqrt(length(costs)),4))
     end
-
-    return stats
 end
 
 
-"""
-Estimate upper bound with Monte Carlo.
-
-# Arguments
-* `model::SPmodel`:
-    the stochastic problem we want to optimize
-* `param::SDDPparameters`:
-    the parameters of the SDDP algorithm
-* `V::Array{PolyhedralFunction}`:
-    the current estimation of Bellman's functions
-* `problems::Array{JuMP.Model}`:
-    Linear model used to approximate each value function
-* `n_simulation::Float64`:
-    Number of scenarios to use to compute Monte-Carlo estimation
-
-# Return
-* `upb::Float64`:
-    estimation of upper bound
-* `costs::Vector{Float64}`:
-    Costs along different trajectories
-"""
-function estimate_upper_bound(model::SPModel, param::SDDPparameters,
-                                V::Vector{PolyhedralFunction},
-                                problem::Vector{JuMP.Model},
-                                n_simulation=1000::Int)
-    aleas = simulate_scenarios(model.noises, n_simulation)
-    return estimate_upper_bound(model, param, aleas, problem)
-end
-function estimate_upper_bound(model::SPModel, param::SDDPparameters,
-                                aleas::Array{Float64, 3},
-                                problem::Vector{JuMP.Model})
-    costs = forward_simulations(model, param, problem, aleas)[1]
-    return upper_bound(costs), costs
-end
-
 
 """
 Build a cut approximating terminal cost with null function
@@ -485,81 +463,3 @@ function add_cuts_to_model!(model::SPModel, t::Int64, problem::JuMP.Model, V::Po
         @constraint(problem, V.betas[i] + dot(lambda, xf) <= alpha)
     end
 end
-
-
-"""
-Exact pruning of all polyhedral functions in input array.
-
-# Arguments
-* `model::SPModel`:
-* `params::SDDPparameters`:
-* `Vector{PolyhedralFunction}`:
-    Polyhedral functions where cuts will be removed
-"""
-function prune_cuts!(model::SPModel, params::SDDPparameters, V::Vector{PolyhedralFunction})
-    for i in 1:length(V)-1
-        V[i] = exact_prune_cuts(model, params, V[i])
-    end
-end
-
-
-"""
-Remove useless cuts in PolyhedralFunction.
-
-# Arguments
-* `model::SPModel`:
-* `params::SDDPparameters`:
-* `V::PolyhedralFunction`:
-    Polyhedral function where cuts will be removed
-
-# Return
-* `PolyhedralFunction`: pruned polyhedral function
-"""
-function exact_prune_cuts(model::SPModel, params::SDDPparameters, V::PolyhedralFunction)
-    ncuts = V.numCuts
-    # Find all active cuts:
-    if ncuts > 1
-        active_cuts = Bool[is_cut_relevant(model, i, V, params.SOLVER) for i=1:ncuts]
-        return PolyhedralFunction(V.betas[active_cuts], V.lambdas[active_cuts, :], sum(active_cuts))
-    else
-        return V
-    end
-end
-
-
-"""
-Test whether the cut number k is relevant to define polyhedral function Vt.
-
-# Arguments
-* `model::SPModel`:
-* `k::Int`:
-    Position of cut to test in PolyhedralFunction object
-* `Vt::PolyhedralFunction`:
-    Object storing all cuts
-* `solver`:
-    Solver to use to solve linear problem
-* `epsilon::Float64`: default is `1e-5`
-    Acceptable tolerance to test cuts relevantness
-
-# Return
-* `Bool`: true if the cut is useful in the definition, false otherwise
-"""
-function is_cut_relevant(model::SPModel, k::Int, Vt::PolyhedralFunction, solver; epsilon=1e-5)
-    m = Model(solver=solver)
-    @variable(m, alpha)
-    @variable(m, model.xlim[i][1] <= x[i=1:model.dimStates] <= model.xlim[i][2])
-
-    for i in 1:Vt.numCuts
-        if i!=k
-            lambda = vec(Vt.lambdas[i, :])
-            @constraint(m, Vt.betas[i] + dot(lambda, x) <= alpha)
-        end
-    end
-
-    λ_k = vec(Vt.lambdas[k, :])
-    @objective(m, Min, alpha - dot(λ_k, x) - Vt.betas[k])
-    solve(m)
-    sol = getobjectivevalue(m)
-    return sol < epsilon
-end
-

src/StochDynamicProgramming.jl

@@ -29,4 +29,6 @@ include("SDDPoptimize.jl")
 include("extensiveFormulation.jl")
 include("SDPoptimize.jl")
 include("compare.jl")
+include("cutpruning.jl")
+include("stoppingtest.jl")
 end