diff --git a/.travis.yml b/.travis.yml index 5321b7a..605b101 100644 --- a/.travis.yml +++ b/.travis.yml @@ -9,6 +9,9 @@ notifications: script: - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi - julia -e 'Pkg.clone(pwd()); Pkg.build("StochDynamicProgramming"); Pkg.test("StochDynamicProgramming"; coverage=true)' +before_script: + - julia -e 'Pkg.clone("https://github.com/JuliaPolyhedra/CutPruners.jl")' + - julia -e 'Pkg.clone("https://github.com/blegat/StochasticDualDynamicProgramming.jl")' after_success: - echo $TRAVIS_JULIA_VERSION - julia -e 'Pkg.add("Coverage"); cd(Pkg.dir("StochDynamicProgramming")); using Coverage; Coveralls.submit(Coveralls.process_folder()); Codecov.submit(Codecov.process_folder())' diff --git a/README.md b/README.md index 43d9bb4..fdcf76b 100644 --- a/README.md +++ b/README.md @@ -28,13 +28,13 @@ It is built upon [JuMP] - Linear dynamics - Linear or convex piecewise linear cost -Extension to non-linear formulation are under development. +Extension to non-linear formulation are under development. Extension to more complex alea dependance are under developpment. ## Why Extensive formulation ? An extensive formulation approach consists in representing the stochastic problem as a deterministic -one with more variable and call a standard deterministic solver. Mainly usable in a linear +one with more variable and call a standard deterministic solver. Mainly usable in a linear setting. Computational complexity is exponential in the number of stages. ## Why Stochastic Dynamic Programming ? @@ -54,10 +54,15 @@ control strategies. ## Installation -Installing StochDynamicProgramming is an easy process. Open Julia and enter +Installing StochDynamicProgramming is an easy process. +Currently, the package depends upon `CutPruners.jl`, which is not +yet registered in Julia's METADATA. To install the package, +open Julia and enter ```julia julia> Pkg.update() +julia> Pkg.clone("https://github.com/JuliaPolyhedra/CutPruners.jl") +julia> Pkg.clone("https://github.com/blegat/StochasticDualDynamicProgramming.jl") julia> Pkg.add("StochDynamicProgramming") ``` diff --git a/REQUIRE b/REQUIRE index f20beb6..6bcd14c 100644 --- a/REQUIRE +++ b/REQUIRE @@ -4,3 +4,5 @@ Distributions ProgressMeter Interpolations Iterators +CutPruners +StochasticDualDynamicProgramming diff --git a/TODO.md b/TODO.md deleted file mode 100644 index 6fb220c..0000000 --- a/TODO.md +++ /dev/null @@ -1,21 +0,0 @@ -# TODO: - - -## SDDP algorithms -- [ ] Add parallelization to forward and backward phase -- [ ] None - -## Examples -- [ ] Merge newsvendor to master -- [ ] Build a utility management example, with hydro-valley and thermal plants -- [ ] Make a benchmark - -## Implementation -- [ ] Fix packaging and importation -- [ ] Ensure that all arrays are read in column-order -- [ ] Improve coverage of unittests - -## Documentation -- [ ] Write examples in `.rst` -- [ ] Write IJulia Notebook -- [ ] Document code in `rst` diff --git a/examples/benchmark.jl b/examples/benchmark.jl deleted file mode 100644 index 03dab8c..0000000 --- a/examples/benchmark.jl +++ /dev/null @@ -1,344 +0,0 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at http://mozilla.org/MPL/2.0/. -############################################################################# -# Benchmark SDDP and DP algorithms upon damsvalley example -# This benchmark includes: -# - comparison of execution time -# - gap between the stochastic solution and the anticipative solution -# -# Usage: -# ``` -# $ julia benchmark.jl {DP|SDDP} -# -# ``` -############################################################################# - -srand(2713) - -using StochDynamicProgramming, JuMP -using Clp - -SOLVER = ClpSolver() -#= const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0, CPX_PARAM_THREADS=4) =# - -const N_STAGES = 10 - -# COST: -const COST = -66*2.7*(1 + .5*(rand(N_STAGES) - .5)) - -# Define dynamic of the dam: -function dynamic(t, x, u, w) - return [x[1] - u[1] - u[3] + w[1], x[2] - u[2] - u[4] + u[1] + u[3]] -end - -# Define cost corresponding to each timestep: -function cost_t(t, x, u, w) - return COST[t] * (u[1] + u[2]) -end - -function final_cost(x) - return 0. -end - - - -"""Solve the problem with a solver, supposing the aleas are known -in advance.""" -function solve_anticipative_problem(model, scenario) - m = Model(solver=SOLVER) - - - @variable(m, model.xlim[1][1] <= x1[1:(N_STAGES)] <= model.xlim[1][2]) - @variable(m, model.xlim[2][1] <= x2[1:(N_STAGES)] <= model.xlim[2][2]) - @variable(m, model.ulim[1][1] <= u1[1:N_STAGES-1] <= model.ulim[1][2]) - @variable(m, model.ulim[2][1] <= u2[1:N_STAGES-1] <= model.ulim[2][2]) - - @objective(m, Min, sum{COST[i]*(u1[i] + u2[i]), i = 1:N_STAGES-1}) - - for i in 1:N_STAGES-1 - @constraint(m, x1[i+1] - x1[i] + u1[i] - scenario[i] == 0) - @constraint(m, x2[i+1] - x2[i] + u2[i] - u1[i] == 0) - end - - @constraint(m, x1[1] == model.initialState[1]) - @constraint(m, x2[1] == model.initialState[2]) - - status = solve(m) - return getobjectivevalue(m) -end - - -"""Build aleas probabilities for each week. - -Return an array with size (N_ALEAS, N_STAGES) - -""" -function build_aleas() - W_MAX = round(Int, .5/7. * 100) - W_MIN = 0 - DW = 1 - # Define aleas' space: - N_ALEAS = Int(round(Int, (W_MAX - W_MIN) / DW + 1)) - ALEAS = linspace(W_MIN, W_MAX, N_ALEAS) - - aleas = zeros(N_ALEAS, N_STAGES) - - # take into account seasonality effects: - unorm_prob = linspace(1, N_ALEAS, N_ALEAS) - proba1 = unorm_prob / sum(unorm_prob) - proba2 = proba1[N_ALEAS:-1:1] - - for t in 1:N_STAGES - aleas[:, t] = (1 - sin(pi*t/N_STAGES)) * proba1 + sin(pi*t/N_STAGES) * proba2 - end - return aleas -end - - -"""Build an admissible scenario for water inflow. - -Parameters: -- n_scenarios (Int64) - Number of scenarios to generate - -- probabilities (Array{Float64, 2}) - Probabilities of occurence of water inflow for each week - -Return an array with size (n_scenarios, N_STAGES) - -""" -function build_scenarios(n_scenarios::Int64, probabilities::Array{Float64}) - scenarios = zeros(n_scenarios, N_STAGES) - - for scen in 1:n_scenarios - for t in 1:N_STAGES - Pcum = cumsum(probabilities[:, t]) - - n_random = rand() - prob = findfirst(x -> x > n_random, Pcum) - scenarios[scen, t] = prob - end - end - return scenarios -end - - -"""Build n_scenarios according to a given distribution. - -Return a Vector{NoiseLaw}""" -function generate_probability_laws(n_scenarios) - aleas = build_scenarios(n_scenarios, build_aleas()) - - laws = Vector{NoiseLaw}(N_STAGES) - - # uniform probabilities: - proba = 1/n_scenarios*ones(n_scenarios) - - for t=1:N_STAGES - laws[t] = NoiseLaw(aleas[:, t], proba) - end - - return laws -end - - -"""Instantiate the problem. - -Return a tuple (SPModel, SDDPparameters) - -""" -function init_problem() - - N_SCENARIOS = 10 - - x0 = [50, 50] - aleas = generate_probability_laws(N_SCENARIOS) - - # Constants: - VOLUME_MAX = 100 - VOLUME_MIN = 0 - CONTROL_MAX = round(Int, .4/7. * VOLUME_MAX) + 1 - CONTROL_MIN = 0 - - - - x_bounds = [(VOLUME_MIN, VOLUME_MAX), (VOLUME_MIN, VOLUME_MAX)] - u_bounds = [(CONTROL_MIN, CONTROL_MAX), (CONTROL_MIN, CONTROL_MAX), (0, Inf), (0, Inf)] - - model = LinearSPModel(N_STAGES, - u_bounds, - x0, - cost_t, - dynamic, - aleas) - - set_state_bounds(model, x_bounds) - - - EPSILON = .05 - MAX_ITER = 20 - solver = SOLVER - params = SDDPparameters(solver, - passnumber=N_SCENARIOS, - gap=EPSILON, - max_iterations=MAX_ITER) - - return model, params -end - - -"""Benchmark SDDP.""" -function benchmark_sddp() - model, params = init_problem() - - # Launch a first start to compile solve_SDDP - params.maxItNumber = 2 - V, pbs = solve_SDDP(model, params, 0) - params.maxItNumber = 20 - # Launch benchmark - println("Launch SDDP ...") - tic() - V, pbs = @time solve_SDDP(model, params, 0) - texec = toq() - println("Time to solve SDDP: ", texec, "s") - - # Test results upon 100 assessment scenarios: - n_assessments = 100 - aleas = simulate_scenarios(model.noises, - n_assessments) - - tic() - costs_sddp, stocks = forward_simulations(model, params, pbs, aleas) - texec = toq() - println("Time to perform simulation: ", texec, "s") - - println("SDDP cost: \t", mean(costs_sddp)) - return stocks, V -end - - -"""Benchmark SDP.""" -function benchmark_sdp() - - N_STAGES::Int64 = 5 - TF = N_STAGES - # Capacity of dams: - VOLUME_MAX::Float64 = 20. - VOLUME_MIN::Float64 = 0. - - # Specify the maximum flow of turbines: - CONTROL_MAX::Float64 = 5. - CONTROL_MIN::Float64 = 0. - - # Some statistics about aleas (water inflow): - W_MAX::Float64 = 5. - W_MIN::Float64 = 0. - DW::Float64 = 1. - - T0 = 1 - - # Define aleas' space: - N_ALEAS = Int(round(Int, (W_MAX - W_MIN) / DW + 1)) - ALEAS = linspace(W_MIN, W_MAX, N_ALEAS); - - N_CONTROLS = 2; - N_STATES = 2; - N_NOISES = 1; - - infoStruct = "HD" - - COST = 66*2.7*(1 + .5*(rand(TF) - .5)); - - # Define dynamic of the dam: - function dynamic(t, x, u, w) - return [x[1] + u[1] + w[1] - u[2], x[2] - u[1]] - end - - # Define cost corresponding to each timestep: - function cost_t(t, x, u, w) - return COST[t] * u[1] - end - - function constraints(t, x, u, w) - return true - end - - function finalCostFunction(x) - return 0. - end - - """Build admissible scenarios for water inflow over the time horizon.""" - function build_scenarios(n_scenarios::Int64) - scenarios = zeros(n_scenarios, TF) - - for scen in 1:n_scenarios - scenarios[scen, :] = (W_MAX-W_MIN)*rand(TF)+W_MIN - end - return scenarios - end - - """Build probability distribution at each timestep based on N scenarios. - Return a Vector{NoiseLaw}""" - function generate_probability_laws(N_STAGES, N_SCENARIOS) - aleas = zeros(N_SCENARIOS, TF, 1) - aleas[:, :, 1] = build_scenarios(N_SCENARIOS) - - laws = Vector{NoiseLaw}(N_STAGES) - - # uniform probabilities: - proba = 1/N_SCENARIOS*ones(N_SCENARIOS) - - for t=1:N_STAGES - aleas_t = reshape(aleas[:, t, :], N_SCENARIOS, 1)' - laws[t] = NoiseLaw(aleas_t, proba) - end - - return laws - end - - N_SCENARIO = 5 - aleas = generate_probability_laws(TF, N_SCENARIO) - - x_bounds = [(VOLUME_MIN, VOLUME_MAX), (VOLUME_MIN, VOLUME_MAX)]; - u_bounds = [(CONTROL_MIN, CONTROL_MAX), (VOLUME_MIN, 10)]; - - x0 = [20., 22.] - - modelSDP = StochDynProgModel(N_STAGES-1, - x_bounds, u_bounds, - x0, cost_t, - finalCostFunction, dynamic, - constraints, aleas); - - stateSteps = [1.,1.]; - controlSteps = [1.,1.]; - monteCarloSize = 10.; - - paramsSDP = StochDynamicProgramming.SDPparameters(modelSDP, stateSteps, - controlSteps, - infoStruct, - "Exact"); - - print("Problem size (T*X*U*W): ") - println(paramsSDP.totalStateSpaceSize*paramsSDP.totalControlSpaceSize) - - n_benchmark = 10 - timing = zeros(n_benchmark) - for n in 1:n_benchmark - tic() - V_sdp = solve_DP(modelSDP, paramsSDP,0); - timing[n] = toq() - end - @show timing - println("SDP execution time: ", mean(timing[2:end]), " s +/-", 3.std(timing[2:end])) - -end - -# SDDP benchmark: -if ARGS[1] == "SDDP" - benchmark_sddp() -elseif ARGS[1] == "DP" - benchmark_sdp() -end diff --git a/examples/dam.jl b/examples/dam.jl index de3549c..6bebe73 100644 --- a/examples/dam.jl +++ b/examples/dam.jl @@ -8,10 +8,9 @@ ############################################################################# -using StochDynamicProgramming, JuMP, Clp - -const SOLVER = ClpSolver() +using StochDynamicProgramming, JuMP +include("solver.jl") const EPSILON = .05 const MAX_ITER = 20 @@ -164,10 +163,10 @@ end function solve_dams(display=0) model, params = init_problem() - V, pbs = solve_SDDP(model, params, display) + sddp = solve_SDDP(model, params, display) aleas = simulate_scenarios(model.noises, params.forwardPassNumber) - costs, stocks = forward_simulations(model, params, pbs, aleas) + costs, stocks = forward_simulations(model, params, sddp.solverinterface, aleas) println("SDDP cost: ", costs) return stocks end diff --git a/examples/damsvalley.jl b/examples/damsvalley.jl index 2d7aaca..d76402e 100644 --- a/examples/damsvalley.jl +++ b/examples/damsvalley.jl @@ -6,7 +6,9 @@ # Set a seed for reproductability: srand(2713) -using StochDynamicProgramming, JuMP, CPLEX +using StochDynamicProgramming, JuMP + +include("solver.jl") ################################################## @@ -116,7 +118,7 @@ function init_problem() set_state_bounds(model, x_bounds) # We need to use CPLEX to solve QP at final stages: - solver = CPLEX.CplexSolver(CPX_PARAM_SIMDISPLAY=0, CPX_PARAM_BARDISPLAY=0) + solver = SOLVER params = SDDPparameters(solver, passnumber=FORWARD_PASS, @@ -128,5 +130,5 @@ end # Solve the problem: model, params = init_problem() -V, pbs = @time solve_SDDP(model, params, 1) +sddp = @time solve_SDDP(model, params, 1) diff --git a/examples/multistock-example.jl b/examples/multistock-example.jl index a08b6e5..c606e3d 100644 --- a/examples/multistock-example.jl +++ b/examples/multistock-example.jl @@ -83,8 +83,9 @@ if run_sddp paramSDDP = SDDPparameters(SOLVER, passnumber=10, max_iterations=MAX_ITER) - V, pbs = solve_SDDP(spmodel, paramSDDP, 2) # display information every 2 iterations - lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, V) + sddp = solve_SDDP(spmodel, paramSDDP, 2, # display information every 2 iterations + stopcrit=StochasticDualDynamicProgramming.IterLimit(MAX_ITER)) + lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions) println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4))) toc(); println(); end @@ -111,7 +112,7 @@ end #srand(1234) # to fix the random seed accross runs if run_sddp && run_sdp && test_simulation scenarios = StochDynamicProgramming.simulate_scenarios(xi_laws,1000) - costsddp, stocks = forward_simulations(spmodel, paramSDDP, pbs, scenarios) + costsddp, stocks = forward_simulations(spmodel, paramSDDP, sddp.solverinterface, scenarios) costsdp, states, controls = sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs) println("Simulated relative gain of sddp over sdp: " *string(round(200*mean(costsdp-costsddp)/abs(mean(costsddp+costsdp)),3))*"%") diff --git a/src/SDDPoptimize.jl b/src/SDDPoptimize.jl index 65928d1..2dbdb4e 100644 --- a/src/SDDPoptimize.jl +++ b/src/SDDPoptimize.jl @@ -1,4 +1,4 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. @@ -9,8 +9,12 @@ ############################################################################# +export solve_SDDP, solve! + """ -Solve spmodel using SDDP algorithm and return lower approximation of Bellman functions. +Solve spmodel using SDDP algorithm and return `SDDPInterface` instance. + +$(SIGNATURES) # Description Alternate forward and backward phases untill the stopping criterion is @@ -27,28 +31,29 @@ fulfilled. `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`: statistics of the algorithm run +`SDDPInterface` """ -function solve_SDDP(model::SPModel, param::SDDPparameters, verbose=0::Int64) - check_SDDPparameters(model,param,verbose) - # initialize value functions: - V, problems = initialize_value_functions(model, param) - (verbose > 0) && println("Initial value function loaded into memory.") +function solve_SDDP(model::SPModel, param::SDDPparameters, verbose=0::Int64; + stopcrit::AbstractStoppingCriterion=IterLimit(20), + prunalgo::AbstractCutPruningAlgo=CutPruners.AvgCutPruningAlgo(-1), + regularization=nothing) # Run SDDP: - sddp_stats = run_SDDP!(model, param, V, problems, verbose) - return V, problems, sddp_stats + sddp = SDDPInterface(model, param, + stopcrit, + prunalgo, + verbose=verbose, + regularization=regularization) + solve!(sddp) + sddp end """ -Solve spmodel using SDDP algorithm and return lower approximation of Bellman functions. +Solve spmodel using SDDP algorithm and return `SDDPInterface` instance. Use hotstart. +$(SIGNATURES) + # Description Alternate forward and backward phases untill the stopping criterion is fulfilled. @@ -66,118 +71,114 @@ fulfilled. `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`: statistics of the algorithm run +* `SDDPInterface` """ -function solve_SDDP(model::SPModel, param::SDDPparameters, V::Vector{PolyhedralFunction}, verbose=0::Int64) - 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 +function solve_SDDP(model::SPModel, param::SDDPparameters, V::Vector{PolyhedralFunction}, verbose=0::Int64; + stopcrit::AbstractStoppingCriterion=IterLimit(20), + prunalgo::AbstractCutPruningAlgo=CutPruners.AvgCutPruningAlgo(-1)) + sddp = SDDPInterface(model, param, + stopcrit, + prunalgo, V, + verbose=verbose) + solve!(sddp) + sddp end -"""Run SDDP iterations. +"""Run SDDP iterations on `sddp::SDDPInterface` instance. -# 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. +$(SIGNATURES) + +# Description +This function modifies `sddp`: +* if `sddp.init` is false, init `sddp` +* run SDDP iterations and update `sddp` till stopping test is fulfilled + +At each iteration, the algorithm runs: +* a forward pass on `sddp` to compute `trajectories` +* a backward pass to update value functions of `sddp` +* a cut pruning to remove outdated cuts in `sddp` +* an estimation of the upper-bound of `sddp` +* an update of the different attributes of `sddp` +* test the stopping criterion -# 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() - - # Initialize cuts container for cuts pruning: - if isa(param.pruning[:type], Union{Type{Territory}, Type{LevelOne}}) - activecuts = [ActiveCutsContainer(model.dimStates) for i in 1:model.stageNumber-1] - else - activecuts = [nothing for i in 1:model.stageNumber-1] - end +function solve!(sddp::SDDPInterface) - (verbose > 0) && println("Initialize cuts") + if ~sddp.init + init!(sddp) + end + model = sddp.spmodel + param = sddp.params + stats = sddp.stats # If computation of upper-bound is needed, a set of scenarios is built # to keep always the same realization for upper bound estimation: - upperbound_scenarios = simulate_scenarios(model.noises, param.in_iter_mc) + upperbound_scenarios = simulate_scenarios(sddp.spmodel.noises, sddp.params.in_iter_mc) upb = [Inf, Inf, Inf] stopping_test::Bool = false # Launch execution of forward and backward passes: - while (~stopping_test) + while !stop(sddp.stopcrit, stats) # Time execution of current pass: tic() #################### # Forward pass : compute stockTrajectories - costs, stockTrajectories, callsolver_forward = forward_pass!(model, param, V, problems) + costs, trajectories = forward_pass!(sddp) #################### # Backward pass : update polyhedral approximation of Bellman functions - callsolver_backward = backward_pass!(model, param, V, problems, stockTrajectories, model.noises) + backward_pass!(sddp, trajectories, model.noises) #################### # Time execution of current pass - lwb = get_bellman_value(model, param, 1, V[1], model.initialState) time_pass = toq() #################### # cut pruning - if param.pruning[:pruning] - prune_cuts!(model, param, V, stockTrajectories, activecuts, stats.niterations+1, verbose) - if (stats.niterations%param.pruning[:period]==0) - problems = hotstart_SDDP(model, param, V) - end - end + prune!(sddp, trajectories) #################### - # In iteration upper bound estimation - upb = in_iteration_upb_estimation(model, param, stats.niterations+1, verbose, - upperbound_scenarios, upb, problems) - - updateSDDPStat!(stats, callsolver_forward+callsolver_backward, lwb, upb, time_pass) - print_current_stats(stats,verbose) + # In iteration lower bound estimation + lwb = lowerbound(sddp) #################### - # Stopping test - stopping_test = test_stopping_criterion(param,stats) + # In iteration upper bound estimation + upb = in_iteration_upb_estimation(model, param, stats.niterations+1, sddp.verbose, + upperbound_scenarios, upb, + sddp.solverinterface) + + updateSDDP!(sddp, lwb, upb, time_pass, trajectories) + checkit(sddp.verbose, sddp.stats.niterations) && println(sddp.stats) end ########## # Estimate final upper bound with param.monteCarloSize simulations: - display_final_solution(model, param, V, problems, stats, verbose) - return stats + finalpass!(sddp) +end + + +"""Init `sddp::SDDPInterface` object.""" +function init!(sddp::SDDPInterface) + random_pass!(sddp) + sddp.init = true + (sddp.verbose > 0) && println("Initialize cuts") end """Display final results once SDDP iterations are finished.""" -function display_final_solution(model::SPModel, param::SDDPparameters, V, - problems, stats::SDDPStat, verbose::Int64) +function finalpass!(sddp::SDDPInterface) + model = sddp.spmodel + param = sddp.params + V = sddp.bellmanfunctions + problems = sddp.solverinterface + stats = sddp.stats + verbose = sddp.verbose + if (verbose>0) && (param.compute_ub >= 0) - lwb = get_bellman_value(model, param, 1, V[1], model.initialState) + lwb = lowerbound(sddp) if (param.compute_ub == 0) || (param.monteCarloSize > 0) (verbose > 0) && println("Compute final upper-bound with ", @@ -189,7 +190,7 @@ function display_final_solution(model::SPModel, param::SDDPparameters, V, σ = stats.upper_bounds_std[end] end - println("\n############################################################") + println("\n", "#"^60) println("SDDP CONVERGENCE") @printf("- Exact lower bound: %.4e [Gap < %.2f%s]\n", lwb, 100*(upb+tol-lwb)/lwb, '%') @@ -197,23 +198,48 @@ function display_final_solution(model::SPModel, param::SDDPparameters, V, @printf("- Upper-bound's s.t.d: %.4e\n", σ) @printf("- Confidence interval (%d%s): [%.4e , %.4e]", 100*(1- 2*(1-param.confidence_level)), '\%',upb-tol, upb+tol) - println("\n############################################################") + println("\n", "#"^60) end end +function updateSDDP!(sddp::SDDPInterface, lwb, upb, time_pass, trajectories) + # Update SDDP stats + updateSDDPStat!(sddp.stats, lwb, upb, time_pass) + + # If specified, reload JuMP model + # this step can be useful if MathProgBase interface takes too much + # room in memory, rendering necessary a call to GC + if checkit(sddp.params.reload, sddp.stats.niterations) + sddp.solverinterface = hotstart_SDDP(sddp.spmodel, + sddp.params, + sddp.bellmanfunctions) + end + + # Update regularization + if !isnull(sddp.regularizer) + update_penalization!(get(sddp.regularizer)) + get(sddp.regularizer).incumbents = trajectories + end +end + """ -Build a cut approximating terminal cost with null function +Build final cost with PolyhedralFunction function `Vt`. + +$(SIGNATURES) # Arguments +* `model::SPModel`: + Model description * `problem::JuMP.Model`: Cut approximating the terminal cost -* `shape`: - If PolyhedralFunction is given, build terminal cost with it - Else, terminal cost is null +* `Vt::PolyhedralFunction`: + Final cost given as a PolyhedralFunction """ -function build_terminal_cost!(model::SPModel, problem::JuMP.Model, Vt::PolyhedralFunction) +function build_terminal_cost!(model::SPModel, + problem::JuMP.Model, + Vt::PolyhedralFunction) # if shape is PolyhedralFunction, build terminal cost with it: alpha = getvariable(problem, :alpha) xf = getvariable(problem, :xf) @@ -230,7 +256,9 @@ end """ -Initialize each linear problem used to approximate value functions +Initialize each linear problem used to approximate value functions + +$(SIGNATURES) # Description This function define the variables and the constraints of each @@ -268,11 +296,11 @@ function build_model(model, param, t) @constraint(m, xf .== model.dynamics(t, x, u, w)) # Add equality and inequality constraints: - if model.equalityConstraints != nothing - @constraint(m, model.equalityConstraints(t, x, u, w) .== 0) + if ~isnull(model.equalityConstraints) + @constraint(m, get(model.equalityConstraints)(t, x, u, w) .== 0) end - if model.inequalityConstraints != nothing - @constraint(m, model.inequalityConstraints(t, x, u, w) .<= 0) + if ~isnull(model.inequalityConstraints) + @constraint(m, get(model.inequalityConstraints)(t, x, u, w) .<= 0) end # Define objective function (could be linear or piecewise linear) @@ -305,6 +333,8 @@ end """ Initialize value functions along a given trajectory +$(SIGNATURES) + # Description This function add the fist cut to each PolyhedralFunction stored in a Array @@ -324,8 +354,7 @@ function initialize_value_functions(model::SPModel, param::SDDPparameters) solverProblems = build_models(model, param) - V = PolyhedralFunction[ - PolyhedralFunction(model.dimStates) for i in 1:model.stageNumber] + V = getemptyvaluefunctions(model) # Build scenarios according to distribution laws: aleas = simulate_scenarios(model.noises, param.forwardPassNumber) @@ -336,19 +365,33 @@ function initialize_value_functions(model::SPModel, build_terminal_cost!(model, solverProblems[end], V[end]) elseif isa(model.finalCost, Function) # In this case, define a trivial value functions for final cost to avoid problem: - V[end] = PolyhedralFunction(zeros(1), zeros(1, model.dimStates), 1) + V[end] = PolyhedralFunction(zeros(1), zeros(1, model.dimStates), 1, UInt64[], 0) model.finalCost(model, solverProblems[end]) end + return V, solverProblems +end +getemptyvaluefunctions(model) = PolyhedralFunction[PolyhedralFunction(model.dimStates) for i in 1:model.stageNumber] + + +""" +Run SDDP iteration with random forward pass. + +$(SIGNATURES) + +# Parameters +* `sddp:SDDPInterface` + SDDP instance +""" +function random_pass!(sddp::SDDPInterface) + model = sddp.spmodel + param = sddp.params stockTrajectories = zeros(model.stageNumber, param.forwardPassNumber, model.dimStates) for i in 1:model.stageNumber, j in 1:param.forwardPassNumber stockTrajectories[i, j, :] = get_random_state(model) end - callsolver = backward_pass!(model, param, V, solverProblems, - stockTrajectories, model.noises) - - return V, solverProblems + backward_pass!(sddp, stockTrajectories, model.noises) end @@ -356,6 +399,8 @@ end Initialize JuMP.Model vector with a previously computed PolyhedralFunction vector. +$(SIGNATURES) + # Arguments * `model::SPModel`: Parametrization of the problem @@ -387,7 +432,9 @@ end """ -Compute value of Bellman function at point xt. Return V_t(xt) +Compute value of Bellman function at point `xt`. Return `V_t(xt)`. + +$(SIGNATURES) # Arguments * `model::SPModel`: @@ -420,10 +467,24 @@ function get_bellman_value(model::SPModel, param::SDDPparameters, return getvalue(alpha) end +""" +Get lower bound of SDDP instance `sddp`. + +$(SIGNATURES) + +""" +function lowerbound(sddp::SDDPInterface) + return get_bellman_value(sddp.spmodel, sddp.params, 1, + sddp.bellmanfunctions[1], + sddp.spmodel.initialState) +end + """ Compute lower-bound of the problem at initial time. +$(SIGNATURES) + # Arguments * `model::SPModel`: Parametrization of the problem @@ -444,6 +505,8 @@ end """ Compute optimal control at point xt and time t. +$(SIGNATURES) + # Arguments * `model::SPModel`: Parametrization of the problem @@ -463,13 +526,15 @@ Compute optimal control at point xt and time t. """ function get_control(model::SPModel, param::SDDPparameters, lpproblem::Vector{JuMP.Model}, t::Int, xt::Vector{Float64}, xi::Vector{Float64}) - return solve_one_step_one_alea(model, param, lpproblem[t], t, xt, xi)[2].optimal_control + return solve_one_step_one_alea(model, param, lpproblem[t], t, xt, xi)[1].uopt end """ Add several cuts to JuMP.Model from a PolyhedralFunction +$(SIGNATURES) + # Arguments * `model::SPModel`: Store the problem definition diff --git a/src/SDPoptimize.jl b/src/SDPoptimize.jl index 584c107..db212ec 100644 --- a/src/SDPoptimize.jl +++ b/src/SDPoptimize.jl @@ -1,11 +1,9 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud, Henri Gerard and -# Tristan Rigaut +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # Stochastic dynamic programming algorithm -# ############################################################################# using ProgressMeter diff --git a/src/SDPutils.jl b/src/SDPutils.jl index 59cb4ad..2a7d01a 100644 --- a/src/SDPutils.jl +++ b/src/SDPutils.jl @@ -1,3 +1,11 @@ +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# Dynamic Programming core functions +############################################################################# + module SDPutils using Interpolations diff --git a/src/StochDynamicProgramming.jl b/src/StochDynamicProgramming.jl index da4a16a..da48c9b 100644 --- a/src/StochDynamicProgramming.jl +++ b/src/StochDynamicProgramming.jl @@ -1,27 +1,35 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # SDDP is an implementation of the Stochastic Dual Dynamic Programming # algorithm for multi-stage stochastic convex optimization problem -# see TODO ############################################################################# + +__precompile__() + module StochDynamicProgramming include("SDPutils.jl") -using MathProgBase, JuMP, Distributions +using MathProgBase, JuMP, Distributions, StochasticDualDynamicProgramming +using DocStringExtensions +using CutPruners export solve_SDDP, NoiseLaw, simulate_scenarios, SDDPparameters, LinearSPModel, set_state_bounds, - extensive_formulation, - PolyhedralFunction, NextStep, forward_simulations, + extensive_formulation, + PolyhedralFunction, forward_simulations, StochDynProgModel, SDPparameters, solve_DP, sdp_forward_simulation, sampling, get_control, get_bellman_value, - benchmark_parameters + benchmark_parameters, SDDPInterface +include("noises.jl") include("objects.jl") +include("params.jl") +include("regularization.jl") +include("interface.jl") include("utils.jl") include("oneStepOneAleaProblem.jl") include("forwardBackwardIterations.jl") diff --git a/src/compare.jl b/src/compare.jl index 2db1e85..a83143e 100644 --- a/src/compare.jl +++ b/src/compare.jl @@ -1,4 +1,4 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. @@ -60,7 +60,7 @@ function benchmark_parameters(model, simulationmemory = m2+m3 g = round(100*(upb-V0)/V0) - push!(solver_calls, sddpstats.ncallsolver) + push!(solver_calls, sddpstats.nsolved) push!(solving_times, t1) push!(solving_mem, m1) push!(gap_sols, g) @@ -72,7 +72,7 @@ function benchmark_parameters(model, print("Simulation time = ",round(simulationtime,4),"\t") print("Simulation memory = ", simulationmemory,"\t") print("Gap < ", g,"% with prob 97.5%\t") - println("number external solver call = ", sddpstats.ncallsolver) + println("number external solver call = ", sddpstats.nsolved) end end return solver_calls, solving_times, solving_mem, gap_sols diff --git a/src/cutpruning.jl b/src/cutpruning.jl index 8bd630a..abcc0ba 100644 --- a/src/cutpruning.jl +++ b/src/cutpruning.jl @@ -1,41 +1,10 @@ ############################################################################# -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# -type ActiveCutsContainer - numCuts::Int - territories::Array{Array} #set of states where cut k is active - nstates::Int - states::Array{Float64, 2} #set of states where cuts are tested -end - - -""" Remove redundant cuts in Polyhedral Value function `V` - -# Arguments -* `V::PolyhedralFunction`: -""" -function remove_cuts!(V::PolyhedralFunction) - Vf = hcat(V.lambdas, V.betas) - Vf = unique(Vf, 1) - return PolyhedralFunction(Vf[:, end], Vf[:, 1:end-1], size(Vf)[1]) -end - - -""" Remove redundant cuts in a vector of Polyhedral Functions `Vts`. - -# Arguments -* `Vts::Vector{PolyhedralFunction}`: -""" -function remove_redundant_cuts!(Vts::Vector{PolyhedralFunction}) - n_functions = length(Vts)-1 - for i in 1:n_functions - Vts[i] = remove_cuts!(Vts[i]) - end -end """ @@ -43,209 +12,25 @@ Exact pruning of all polyhedral functions in input array. # Arguments * `model::SPModel`: -* `param::SDDPparameters`: -* `Vector{PolyhedralFunction}`: - Polyhedral functions where cuts will be removed * `trajectories::Array{Float64, 3}` Previous trajectories -* `territory::Array{ActiveCutsContainer}` - Container storing the territory (i.e. the set of tested states where - a given cut is active) for each cuts -* `it::Int64`: - current iteration number -* `verbose::Int64` """ -function prune_cuts!(model::SPModel, - param::SDDPparameters, - V::Vector{PolyhedralFunction}, - trajectories::Array{Float64, 3}, - territory::Union{Array{Void}, Array{ActiveCutsContainer}}, - it::Int64, - verbose::Int64) +function prune!(sddp::SDDPInterface, + trajectories::Array{Float64, 3}, + ) # Basic pruning: remove redundant cuts - remove_redundant_cuts!(V) - - # If pruning is performed with territory heuristic, update territory - # at given iteration: - if isa(param.pruning[:type], Union{Type{Territory}, Type{LevelOne}}) - for t in 1:model.stageNumber-1 - states = reshape(trajectories[t, :, :], param.forwardPassNumber, model.dimStates) - find_level1_cuts!(territory[t], V[t], states) - end - end - - # If specified to prune cuts at this iteration, do it: - if param.pruning[:pruning] && (it%param.pruning[:period]==0) - # initial number of cuts: - ncuts_initial = get_total_number_cuts(V) - (verbose > 0) && print("Prune cuts ...") - - for i in 1:length(V)-1 - V[i] = pcuts!(param.pruning[:type], model, param, V[i], territory[i]) - end - - # final number of cuts: - ncuts_final = get_total_number_cuts(V) - - (verbose > 0) && @printf(" New cuts/Old cuts ratio: %.3f \n", ncuts_final/ncuts_initial) - end -end - -pcuts!(::Type{LevelOne}, model, param, V, territory) = level1_cuts_pruning!(model, param, V, territory) -pcuts!(::Type{ExactPruning}, model, param, V, territory) = exact_cuts_pruning(model, param, V, territory) -pcuts!(::Type{Territory}, model, param, V, territory) = exact_cuts_pruning_accelerated!(model, param, V, territory) - - -""" -Remove useless cuts in PolyhedralFunction. - -# Arguments -* `model::SPModel`: -* `param::SDDPparameters`: -* `V::PolyhedralFunction`: - Polyhedral function where cuts will be removed - -# Return -* `PolyhedralFunction`: pruned polyhedral function -""" -function exact_cuts_pruning(model::SPModel, param::SDDPparameters, V::PolyhedralFunction, territory) - ncuts = V.numCuts - # Find all active cuts: - if ncuts > 1 - active_cuts = Bool[is_cut_relevant(model, i, V, param.SOLVER)[1] 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), getvalue(x) -end - - -######################################## -# Territory algorithm -######################################## - -ActiveCutsContainer(ndim) = ActiveCutsContainer(0, [], 0, Array{Float64}(0, ndim)) - + #= for t in 1:sddp.spmodel.stageNumber-1 =# + #= b, A = fetchnewcuts!(sddp.bellmanfunctions[t]) =# + #= nnew = length(b) =# + #= if nnew > 0 =# + #= mycut = Bool[true for _ in 1:length(b)] =# + #= CutPruners.addcuts!(sddp.pruner[t], A, b, mycut) =# + #= end =# + #= end =# -"""Update territories (i.e. the set of tested states where - a given cut is active) with cuts previously computed during backward pass. - -# Arguments -* `cutscontainer::ActiveCutsContainer`: -* `Vt::PolyhedralFunction`: - Object storing all cuts -* `states`: - Object storing all visited states -""" -function find_level1_cuts!(cutscontainer::ActiveCutsContainer, V::PolyhedralFunction, states) - nc = V.numCuts - # get number of new positions to analyse: - nx = size(states, 1) - nt = nc - cutscontainer.numCuts - - for i in 1:nt - add_cut!(cutscontainer) - update_territory!(cutscontainer, V, nc - nt + i) - end - - # ensure that territory has the same number of cuts as V! - assert(cutscontainer.numCuts == V.numCuts) - - for i in 1:nx - x = collect(states[i, :]) - add_state!(cutscontainer, V, x) - end end -"""Update territories (i.e. the set of tested states where - a given cut is active) considering new cut given by index `indcut`. - -# Arguments -* `cutscontainer::ActiveCutsContainer`: -* `V::PolyhedralFunction`: - Object storing all cuts -* `indcut::Int64`: - new cut index -""" -function update_territory!(cutscontainer::ActiveCutsContainer, V::PolyhedralFunction, indcut::Int64) - for k in 1:cutscontainer.numCuts - if k == indcut - continue - end - todelete = [] - for (num, (ix, cost)) in enumerate(cutscontainer.territories[k]) - x = collect(cutscontainer.states[ix, :]) - - costnewcut = cutvalue(V, indcut, x) - - if costnewcut > cost - push!(todelete, num) - push!(cutscontainer.territories[indcut], (ix, costnewcut)) - end - end - deleteat!(cutscontainer.territories[k], todelete) - end -end - - -"""Add cut to `ActiveCutsContainer`.""" -function add_cut!(cutscontainer::ActiveCutsContainer) - push!(cutscontainer.territories, []) - cutscontainer.numCuts += 1 -end - - -"""Add a new state to test and accordingly update territories of each cut.""" -function add_state!(cutscontainer::ActiveCutsContainer, V::PolyhedralFunction, x::Array{Float64}) - # Get cut which is active at point `x`: - bcost, bcuts = optimalcut(x, V) - - # Add `x` to the territory of cut `bcuts`: - cutscontainer.nstates += 1 - push!(cutscontainer.territories[bcuts], (cutscontainer.nstates, bcost)) - - # Add `x` to the list of visited state: - cutscontainer.states = vcat(cutscontainer.states, x') -end - """Remove cuts in PolyhedralFunction that are inactive on all visited states. # Arguments @@ -257,14 +42,12 @@ end the new PolyhedralFunction """ function level1_cuts_pruning!(model::SPModel, param::SDDPparameters, - V::PolyhedralFunction, cutscontainer::ActiveCutsContainer) - assert(cutscontainer.numCuts == V.numCuts) + V::PolyhedralFunction, cutscontainer::CutPruners.LevelOneCutPruner) nstates = [length(terr) for terr in cutscontainer.territories] active_cuts = nstates .> 0 cutscontainer.territories = cutscontainer.territories[active_cuts] - cutscontainer.numCuts = sum(active_cuts) return PolyhedralFunction(V.betas[active_cuts], V.lambdas[active_cuts, :], sum(active_cuts)) @@ -288,7 +71,7 @@ then test remaining cuts. """ function exact_cuts_pruning_accelerated!(model::SPModel, param::SDDPparameters, V::PolyhedralFunction, - cutscontainer::ActiveCutsContainer) + cutscontainer::CutPruners.LevelOneCutPruner) assert(cutscontainer.numCuts == V.numCuts) solver = param.SOLVER @@ -319,62 +102,7 @@ function exact_cuts_pruning_accelerated!(model::SPModel, param::SDDPparameters, end -"""Find active cut at point `xf`. - -# Arguments -* `xf::Vector{Float64}`: -* `V::PolyhedralFunction`: - Object storing all cuts - -# Return -`bestcost::Float64` - Value of maximal cut at point `xf` -`bestcut::Int64` - Index of maximal cut at point `xf` -""" -function optimalcut(xf::Vector{Float64}, V::PolyhedralFunction) - bestcost = -Inf::Float64 - bestcut = -1 - nstates = size(V.lambdas, 2) - ncuts = size(V.lambdas, 1) - - @inbounds for i in 1:ncuts - cost = V.betas[i] - for j in 1:nstates - cost += V.lambdas[i, j]*xf[j] - end - if cost > bestcost - bestcost = cost - bestcut = i - end - end - return bestcost, bestcut -end - - -""" -Get value of cut with index `indc` at point `x`. - -# Arguments -- `V::PolyhedralFunction` - Approximation of the value function as linear cuts -- `indc::Int64` - Index of cut to consider -- `x::Array{Float64}` - Coordinates of state - -# Return -`cost::Float64` - Value of cut `indc` at point `x` -""" -function cutvalue(V::PolyhedralFunction, indc::Int, x::Array{Float64}) - cost = V.betas[indc] - for j in 1:size(V.lambdas, 2) - cost += V.lambdas[indc, j]*x[j] - end - return cost -end - # Return total number of cuts in PolyhedralFunction array: -get_total_number_cuts(V::Array{PolyhedralFunction}) = sum([v.numCuts for v in V]) +ncuts(V::PolyhedralFunction) = length(V.betas) +ncuts(V::Array{PolyhedralFunction}) = sum([ncuts(v) for v in V]) diff --git a/src/extensiveFormulation.jl b/src/extensiveFormulation.jl index abd476a..af8beb2 100644 --- a/src/extensiveFormulation.jl +++ b/src/extensiveFormulation.jl @@ -1,4 +1,4 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. diff --git a/src/forwardBackwardIterations.jl b/src/forwardBackwardIterations.jl index 1883628..12f957c 100644 --- a/src/forwardBackwardIterations.jl +++ b/src/forwardBackwardIterations.jl @@ -7,19 +7,17 @@ ############################################################################# """ -Make a forward pass of the algorithm +Run a forward pass of the algorithm with `sddp` object + +$(SIGNATURES) # Description Simulate scenarios of noise and compute optimal trajectories on those scenarios, with associated costs. # Arguments -* `model::SPmodel`: the stochastic problem we want to optimize -* `param::SDDPparameters`: the parameters of the SDDP algorithm -* `V::Vector{PolyhedralFunction}`: - Linear model used to approximate each value function -* `problems::Vector{JuMP.Model}`: - Current linear problems +* `sddp::SDDPInterface`: + SDDP interface object # Returns * `costs::Array{float,1}`: @@ -27,33 +25,37 @@ scenarios, with associated costs. * `stockTrajectories::Array{float}`: the simulated stock trajectories. stocks(t,k,:) is the stock for scenario k at time t. -* `callsolver_forward::Int64`: - number of call to solver + """ -function forward_pass!(model::SPModel, - param::SDDPparameters, - V::Vector{PolyhedralFunction}, - problems::Vector{JuMP.Model}) +function forward_pass!(sddp::SDDPInterface) + model = sddp.spmodel + param = sddp.params + solverProblems = sddp.solverinterface + V = sddp.bellmanfunctions + problems = sddp.solverinterface # Draw a set of scenarios according to the probability # law specified in model.noises: noise_scenarios = simulate_scenarios(model.noises, param.forwardPassNumber) # If acceleration is ON, need to build a new array of problem to # avoid side effect: - problems_fp = (param.IS_ACCELERATED)? hotstart_SDDP(model, param, V):problems - costs, stockTrajectories,_,callsolver_forward = forward_simulations(model, + problems_fp = isregularized(sddp) ? hotstart_SDDP(model, param, V) : problems + costs, stockTrajectories,_,callsolver_forward, tocfw = forward_simulations(model, param, problems_fp, noise_scenarios, - acceleration=param.IS_ACCELERATED) + regularizer=sddp.regularizer) - model.refTrajectories = stockTrajectories - return costs, stockTrajectories, callsolver_forward + sddp.stats.nsolved += callsolver_forward + sddp.stats.solverexectime_fw = vcat(sddp.stats.solverexectime_fw, tocfw) + return costs, stockTrajectories end """ -Make a forward pass of the algorithm +Simulate a forward pass of the algorithm + +$(SIGNATURES) # Description Simulate a scenario of noise and compute an optimal trajectory on this @@ -86,9 +88,10 @@ function forward_simulations(model::SPModel, param::SDDPparameters, solverProblems::Vector{JuMP.Model}, xi::Array{Float64}; - acceleration=false) + regularizer=Nullable{SDDPRegularization}()) callsolver::Int = 0 + solvertime = Float64[] T = model.stageNumber nb_forward = size(xi)[2] @@ -123,26 +126,27 @@ function forward_simulations(model::SPModel, callsolver += 1 # Solve optimization problem corresponding to current position: - if acceleration && ~isa(model.refTrajectories, Void) - xp = collect(model.refTrajectories[t+1, k, :]) - status, nextstep = solve_one_step_one_alea(model, param, - solverProblems[t], t, state_t, alea_t, xp) + if !isnull(regularizer) && !isa(get(regularizer).incumbents, Void) + reg = get(regularizer) + xp = getincumbent(reg, t, k) + sol, ts = regularize(model, param, reg, + solverProblems[t], t, state_t, alea_t, xp) else - status, nextstep = solve_one_step_one_alea(model, param, - solverProblems[t], t, state_t, alea_t) + sol, ts = solve_one_step_one_alea(model, param, + solverProblems[t], t, state_t, alea_t) end + push!(solvertime, ts) # Check if the problem is effectively solved: - if status + if sol.status # Get the next position: - stockTrajectories[t+1, k, :] = nextstep.next_state + stockTrajectories[t+1, k, :] = sol.xf # the optimal control just computed: - opt_control = nextstep.optimal_control - controls[t, k, :] = opt_control + controls[t, k, :] = sol.uopt # and the current cost: - costs[k] += nextstep.cost - nextstep.cost_to_go + costs[k] += sol.objval - sol.θ if t==T-1 - costs[k] += nextstep.cost_to_go + costs[k] += sol.θ end else # if problem is not properly solved, next position if equal @@ -153,7 +157,7 @@ function forward_simulations(model::SPModel, end end end - return costs, stockTrajectories, controls, callsolver + return costs, stockTrajectories, controls, callsolver, solvertime end @@ -161,6 +165,8 @@ end """ Add to polyhedral function a cut with shape Vt >= beta + +$(SIGNATURES) + # Arguments * `model::SPModel`: Store the problem definition * `t::Int64`: Current time @@ -172,13 +178,18 @@ Add to polyhedral function a cut with shape Vt >= beta + subgradient of the cut to add """ function add_cut!(model::SPModel, - t::Int64, Vt::PolyhedralFunction, - beta::Float64, lambda::Vector{Float64}) - Vt.lambdas = vcat(Vt.lambdas, reshape(lambda, 1, model.dimStates)) + t::Int64, Vt::PolyhedralFunction, + beta::Float64, lambda::Vector{Float64}) + Vt.lambdas = vcat(Vt.lambdas, lambda') Vt.betas = vcat(Vt.betas, beta) + Vt.hashcuts = vcat(Vt.hashcuts, hash(lambda)) Vt.numCuts += 1 + Vt.newcuts += 1 end +function isinside(Vt::PolyhedralFunction, lambda::Vector{Float64}) + hash(lambda) in Vt.hashcuts +end """ Add a cut to the JuMP linear problem. @@ -198,15 +209,15 @@ Add a cut to the JuMP linear problem. function add_cut_to_model!(model::SPModel, problem::JuMP.Model, t::Int64, beta::Float64, lambda::Vector{Float64}) alpha = getvariable(problem, :alpha) - x = getvariable(problem, :x) - u = getvariable(problem, :u) - w = getvariable(problem, :w) - @constraint(problem, beta + dot(lambda, model.dynamics(t, x, u, w)) <= alpha) + xf = getvariable(problem, :xf) + @constraint(problem, beta + dot(lambda, xf) <= alpha) end """ -Make a backward pass of the algorithm +Run a SDDP backward pass on `sddp`. + +$(SIGNATURES) # Description For t:T-1 -> 0, compute a valid cut of the Bellman function @@ -214,28 +225,25 @@ Vt at the state given by stockTrajectories and add them to the current estimation of Vt. # 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 -* `solverProblems::Array{JuMP.Model}`: - Linear model used to approximate each value function +* `sddp::SDDPInterface`: + SDDP instance * `stockTrajectories::Array{Float64,3}`: stockTrajectories[t,k,:] is the vector of stock where the cut is computed for scenario k and time t. * `law::Array{NoiseLaw}`: Conditionnal distributions of perturbation, for each timestep """ -function backward_pass!(model::SPModel, - param::SDDPparameters, - V::Vector{PolyhedralFunction}, - solverProblems::Vector{JuMP.Model}, +function backward_pass!(sddp::SDDPInterface, stockTrajectories::Array{Float64, 3}, law) + model = sddp.spmodel + param = sddp.params + solverProblems = sddp.solverinterface + V = sddp.bellmanfunctions + callsolver::Int = 0 + solvertime = Float64[] T = model.stageNumber nb_forward = size(stockTrajectories)[2] @@ -250,7 +258,7 @@ function backward_pass!(model::SPModel, subgradient_array = zeros(Float64, model.dimStates, law[t].supportSize) # We collect current state: - state_t = collect(stockTrajectories[t, k, :]) + state_t = stockTrajectories[t, k, :] # We will store probabilities in a temporary array. # It is initialized at 0. If all problem are infeasible for # current timestep, then proba remains equal to 0 and not cut is added. @@ -265,16 +273,17 @@ function backward_pass!(model::SPModel, callsolver += 1 # We solve LP problem with current noise and position: - solved, nextstep = solve_one_step_one_alea(model, param, - solverProblems[t], - t, state_t, alea_t, - relaxation=model.IS_SMIP) + sol, ts = solve_one_step_one_alea(model, param, + solverProblems[t], + t, state_t, alea_t, + relaxation=model.IS_SMIP) + push!(solvertime, ts) - if solved + if sol.status # We catch the subgradient λ: - subgradient_array[:, w] = nextstep.sub_gradient + subgradient_array[:, w] = sol.ρe # and the current cost: - costs[w] = nextstep.cost + costs[w] = sol.objval # and as problem is solved we store current proba in array: proba[w] = law[t].proba[w] end @@ -293,13 +302,19 @@ function backward_pass!(model::SPModel, beta = costs_npass - dot(subgradient, state_t) # Add cut to polyhedral function and JuMP model: - add_cut!(model, t, V[t], beta, subgradient) - if t > 1 - add_cut_to_model!(model, solverProblems[t-1], t, beta, subgradient) + if ~isinside(V[t], subgradient) + sddp.stats.nocuts += 1 + add_cut!(model, t, V[t], beta, subgradient) + if t > 1 + add_cut_to_model!(model, solverProblems[t-1], t, beta, subgradient) + end end end end end - return callsolver + + # update stats + sddp.stats.nsolved += callsolver + sddp.stats.solverexectime_bw = vcat(sddp.stats.solverexectime_bw, solvertime) end diff --git a/src/interface.jl b/src/interface.jl new file mode 100644 index 0000000..4b10788 --- /dev/null +++ b/src/interface.jl @@ -0,0 +1,74 @@ +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# SDDP interface +############################################################################# + +type SDDPInterface + init::Bool + # Stochastic model to solve + spmodel::SPModel + # SDDP parameters + params::SDDPparameters + # statistics + stats::SDDPStat + # stopping criterion + stopcrit::AbstractStoppingCriterion + # cut pruner: + pruner::Vector{CutPruners.AbstractCutPruner} + # regularization scheme: + regularizer::Nullable{AbstractRegularization} + + # solution + bellmanfunctions::Vector{PolyhedralFunction} + solverinterface::Vector{JuMP.Model} + + verbose::Int + + # Init SDDP interface + function SDDPInterface(model::SPModel, # SP Model + param::SDDPparameters,# parameters + stopcrit::AbstractStoppingCriterion, + prunalgo::AbstractCutPruningAlgo; + regularization=nothing, + verbose::Int=1) + check_SDDPparameters(model, param, verbose) + # initialize value functions: + V, problems = initialize_value_functions(model, param) + (verbose > 0) && println("SDDP Interface initialized") + + pruner = initpruner(prunalgo, model.stageNumber, model.dimStates) + #Initialization of stats + stats = SDDPStat() + return new(false, model, param, stats, stopcrit, pruner, regularization, V, + problems, verbose) + end + function SDDPInterface(model::SPModel, + params::SDDPparameters, + stopcrit::AbstractStoppingCriterion, + prunalgo::AbstractCutPruningAlgo, + V::Vector{PolyhedralFunction}; + regularization=nothing, + verbose::Int=1) + check_SDDPparameters(model, params, verbose) + # First step: process value functions if hotstart is called + problems = hotstart_SDDP(model, params, V) + pruner = initpruner(prunalgo, model.stageNumber, model.dimStates) + + stats = SDDPStat() + return new(false, model, params, stats, stopcrit, pruner, regularization, + V, problems, verbose) + end +end + + +function initpruner(algo, nstages, ndim) + # Initialize cuts container for cuts pruning: + return [CutPruners.CutPruner{ndim, Float64}(algo, :Max) for i in 1:nstages-1] +end + + +isregularized(sddp::SDDPInterface) = !isnull(sddp.regularizer) + diff --git a/src/noises.jl b/src/noises.jl index 7a17dbf..b1a22c6 100644 --- a/src/noises.jl +++ b/src/noises.jl @@ -1,4 +1,4 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. diff --git a/src/objects.jl b/src/objects.jl index 8f82f9f..dd1a653 100644 --- a/src/objects.jl +++ b/src/objects.jl @@ -1,4 +1,4 @@ -# Copyright 2014, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. @@ -6,7 +6,6 @@ # Define all types used in this module. ############################################################################# -include("noises.jl") abstract SPModel @@ -17,10 +16,19 @@ type PolyhedralFunction lambdas::Array{Float64,2} #lambdas[k,:] is the subgradient # number of cuts: numCuts::Int64 + hashcuts::Vector{UInt64} + newcuts::Int end -PolyhedralFunction(ndim) = PolyhedralFunction([], Array{Float64}(0, ndim), 0) +PolyhedralFunction(ndim::Int) = PolyhedralFunction(Float64[], Array{Float64}(0, ndim), 0, UInt64[], 0) +PolyhedralFunction(beta, lambda) = PolyhedralFunction(beta, lambda, length(beta), UInt64[], 0) +function fetchnewcuts!(V::PolyhedralFunction) + β = V.betas[end-V.newcuts+1:end] + λ = V.lambdas[end-V.newcuts+1:end, :] + V.newcuts = 0 + return β, λ +end type LinearSPModel <: SPModel # problem dimension @@ -43,10 +51,8 @@ type LinearSPModel <: SPModel finalCost::Union{Function, PolyhedralFunction} controlCat::Vector{Symbol} - equalityConstraints::Union{Void, Function} - inequalityConstraints::Union{Void, Function} - - refTrajectories::Union{Void, Array{Float64, 3}} + equalityConstraints::Nullable{Function} + inequalityConstraints::Nullable{Function} IS_SMIP::Bool @@ -70,7 +76,7 @@ type LinearSPModel <: SPModel if isa(Vfinal, Function) || isa(Vfinal, PolyhedralFunction) Vf = Vfinal else - Vf = PolyhedralFunction(zeros(1), zeros(1, dimStates), 1) + Vf = PolyhedralFunction(zeros(1), zeros(1, dimStates), 1, UInt64[], 0) end isbu = isa(control_cat, Vector{Symbol})? control_cat: [:Cont for i in 1:dimStates] @@ -79,7 +85,7 @@ type LinearSPModel <: SPModel xbounds = [(-Inf, Inf) for i=1:dimStates] return new(nstage, dimControls, dimStates, dimNoises, xbounds, ubounds, - x0, cost, dynamic, aleas, Vf, isbu, eqconstr, ineqconstr, nothing, is_smip) + x0, cost, dynamic, aleas, Vf, isbu, eqconstr, ineqconstr, is_smip) end end @@ -140,89 +146,6 @@ type StochDynProgModel <: SPModel end -# Define alias for cuts pruning algorithm: -typealias LevelOne Val{:LevelOne} -typealias ExactPruning Val{:Exact} -typealias Territory Val{:Exact_Plus} -typealias NoPruning Val{:none} - -type SDDPparameters - # Solver used to solve LP - SOLVER::MathProgBase.AbstractMathProgSolver - # Solver used to solve MILP (default is nothing): - MIPSOLVER::Union{Void, MathProgBase.AbstractMathProgSolver} - # number of scenarios in the forward pass - forwardPassNumber::Int64 - # Admissible gap between lower and upper-bound: - gap::Float64 - # tolerance upon confidence interval: - confidence_level::Float64 - # Maximum iterations of the SDDP algorithms: - maxItNumber::Int64 - # Define the pruning method - pruning::Dict{Symbol, Any} - # Estimate upper-bound every %% iterations: - compute_ub::Int64 - # Number of MonteCarlo simulation to perform to estimate upper-bound: - monteCarloSize::Int64 - # Number of MonteCarlo simulation to estimate the upper bound during one iteration - in_iter_mc::Int64 - # specify whether SDDP is accelerated - IS_ACCELERATED::Bool - # ... and acceleration parameters: - acceleration::Dict{Symbol, Float64} - - function SDDPparameters(solver; passnumber=10, gap=0., confidence=.975, - max_iterations=20, prune_cuts=0, - pruning_algo="none", - compute_ub=-1, montecarlo_final=1000, montecarlo_in_iter=100, - mipsolver=nothing, - rho0=0., alpha=1.) - - if ~(pruning_algo ∈ ["none", "exact+", "level1", "exact"]) - throw(ArgumentError("`pruning_algo` must be `none`, `level1`, `exact` or `exact+`")) - end - is_acc = (rho0 > 0.) - accparams = is_acc? Dict(:ρ0=>rho0, :alpha=>alpha, :rho=>rho0): Dict() - - pruning_algo = (prune_cuts>0)? pruning_algo: "none" - prune_cuts = (pruning_algo != "none")? prune_cuts: 0 - - corresp = Dict("none"=>NoPruning, - "level1"=>LevelOne, - "exact+"=>Territory, - "exact"=>ExactPruning) - prune_cuts = Dict(:pruning=>prune_cuts>0, - :period=>prune_cuts, - :type=>corresp[pruning_algo]) - return new(solver, mipsolver, passnumber, gap, confidence, - max_iterations, prune_cuts, compute_ub, - montecarlo_final, montecarlo_in_iter, is_acc, accparams) - end -end - -""" -Test compatibility of parameters. - -# Arguments -* `model::SPModel`: - Parametrization of the problem -* `param::SDDPparameters`: - Parameters of SDDP -* `verbose:Int64`: - -# Return -`Bool` -""" -function check_SDDPparameters(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 - (model.IS_SMIP && param.IS_ACCELERATED) && error("Acceleration of SMIP not supported") - (verbose > 0) && (model.IS_SMIP) && println("SMIP SDDP") - (verbose > 0) && (param.IS_ACCELERATED) && println("Acceleration: ON") -end - type SDPparameters stateSteps @@ -262,15 +185,12 @@ type SDPparameters end -function set_max_iterations(param::SDDPparameters, n_iter::Int) - param.maxItNumber = n_iter -end # Define an object to store evolution of solution # along iterations: -type SDDPStat +type SDDPStat <: AbstractSDDPStats # Number of iterations: - niterations::Int64 + niterations::Int # evolution of lower bound: lower_bounds::Vector{Float64} # evolution of upper bound: @@ -281,11 +201,26 @@ type SDDPStat upper_bounds_tol::Vector{Float64} # evolution of execution time: exectime::Vector{Float64} + # time used to solve each LP: + solverexectime_fw::Vector{Float64} + solverexectime_bw::Vector{Float64} # number of calls to solver: - ncallsolver::Int64 + nsolved::Int + # number of optimality cuts + nocuts::Int + npaths::Int + # current lower bound + lowerbound::Float64 + # current lower bound + upperbound::Float64 + # upper-bound std: + σ_UB::Float64 + # total time + time::Float64 end -SDDPStat() = SDDPStat(0, [], [], [], [], [], 0) + +SDDPStat() = SDDPStat(0, [], [], [], [], [], [], [], 0, 0, 0, 0., 0., 0., 0.) """ Update the SDDPStat object with the results of current iterations. @@ -302,25 +237,34 @@ Update the SDDPStat object with the results of current iterations. * `time` """ function updateSDDPStat!(stats::SDDPStat, - callsolver_at_it::Int64, lwb::Float64, upb::Vector{Float64}, time) - stats.ncallsolver += callsolver_at_it stats.niterations += 1 push!(stats.lower_bounds, lwb) push!(stats.upper_bounds, upb[1]) push!(stats.upper_bounds_tol, upb[3]) push!(stats.upper_bounds_std, upb[2]) push!(stats.exectime, time) + stats.lowerbound = lwb + stats.upperbound = upb[1] + stats.σ_UB = upb[2] + stats.time += time end -type NextStep - next_state::Array{Float64, 1} - optimal_control::Array{Float64, 1} - sub_gradient - cost::Float64 - cost_to_go::Float64 +type NLDSSolution + # solver status: + status::Bool + # cost: + objval::Float64 + # next position: + xf::Array{Float64, 1} + # optimal control: + uopt::Array{Float64, 1} + # Subgradient: + ρe::Array{Float64, 1} + # cost-to-go: + θ::Float64 end diff --git a/src/oneStepOneAleaProblem.jl b/src/oneStepOneAleaProblem.jl index c2df7b9..91aeb01 100644 --- a/src/oneStepOneAleaProblem.jl +++ b/src/oneStepOneAleaProblem.jl @@ -1,4 +1,4 @@ -# Copyright 2014, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. @@ -73,35 +73,36 @@ function solve_one_step_one_alea(model, solved = (status == :Optimal) end + solvetime = try getsolvetime(m) catch 0 end + if solved optimalControl = getvalue(u) # Return object storing results: - λ = (~model.IS_SMIP || relaxation)? Float64[getdual(m.ext[:cons][i]) for i in 1:model.dimStates]:nothing - - result = NextStep( + result = NLDSSolution( + solved, + getobjectivevalue(m), model.dynamics(t, xt, optimalControl, xi), optimalControl, - λ, - getobjectivevalue(m), + getdual(m.ext[:cons]), getvalue(alpha)) else # If no solution is found, then return nothing result = nothing end - return solved, result + return result, solvetime end # Solve local problem with a quadratic penalization: -function solve_one_step_one_alea(model, param, m::JuMP.Model, t::Int64, - xt::Vector{Float64}, xi::Vector{Float64}, xp::Vector{Float64}) +function regularize(model, param, + regularizer::AbstractRegularization, + m::JuMP.Model, + t::Int64, + xt::Vector{Float64}, xi::Vector{Float64}, xp::Vector{Float64}) # store current objective: pobj = m.obj xf = getvariable(m, :xf) - # copy JuMP model to avoid side effect: - rho = param.acceleration[:rho] - # build quadratic penalty term: - qexp = QuadExpr(rho*dot(xf - xp, xf - xp)) + qexp = getpenaltyexpr(regularizer, xf, xp) # and update model objective: @objective(m, :Min, m.obj + qexp) res = solve_one_step_one_alea(model, param, m, t, xt, xi) @@ -119,7 +120,7 @@ end """Solve original MILP problem.""" function solve_mip!(m, param) - setsolver(m, param.MIPSOLVER) + setsolver(m, get(param.MIPSOLVER)) status = solve(m, relaxation=false) return status == :Optimal end diff --git a/src/params.jl b/src/params.jl new file mode 100644 index 0000000..917c377 --- /dev/null +++ b/src/params.jl @@ -0,0 +1,60 @@ +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# Definition of SDDP parameters +############################################################################# + +type SDDPparameters + # Solver used to solve LP + SOLVER::MathProgBase.AbstractMathProgSolver + # Solver used to solve MILP (default is nothing): + MIPSOLVER::Nullable{MathProgBase.AbstractMathProgSolver} + # number of scenarios in the forward pass + forwardPassNumber::Int64 + # tolerance upon confidence interval: + confidence_level::Float64 + # Estimate upper-bound every %% iterations: + compute_ub::Int64 + # Number of MonteCarlo simulation to perform to estimate upper-bound: + monteCarloSize::Int64 + # Number of MonteCarlo simulation to estimate the upper bound during one iteration + in_iter_mc::Int64 + # Refresh JuMP Model: + reload::Int + + function SDDPparameters(solver; passnumber=10, gap=0., confidence=.975, + max_iterations=20, prune_cuts=0, + pruning_algo="none", + compute_ub=-1, montecarlo_final=1000, montecarlo_in_iter=100, + mipsolver=nothing, + rho0=0., alpha=1., reload=-1) + + return new(solver, mipsolver, passnumber, confidence, + compute_ub, montecarlo_final, montecarlo_in_iter, reload) + end +end + + +""" +Test compatibility of parameters. + +# Arguments +* `model::SPModel`: + Parametrization of the problem +* `param::SDDPparameters`: + Parameters of SDDP +* `verbose:Int64`: + +# Return +`Bool` +""" +function check_SDDPparameters(model::SPModel, param::SDDPparameters, verbose=0::Int64) + if model.IS_SMIP && isnull(param.MIPSOLVER) + error("MIP solver is not defined. Please set `param.MIPSOLVER`") + end + + (verbose > 0) && (model.IS_SMIP) && println("SMIP SDDP") +end + diff --git a/src/regularization.jl b/src/regularization.jl new file mode 100644 index 0000000..4365895 --- /dev/null +++ b/src/regularization.jl @@ -0,0 +1,39 @@ +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# SDDP regularization +############################################################################# + + +export SDDPRegularization + +abstract AbstractRegularization + + +type SDDPRegularization <: AbstractRegularization + ρ::Float64 + alpha::Float64 + incumbents + decay::Float64 + function SDDPRegularization(rho0::Float64, alpha::Float64; decay=1.) + return new(rho0, alpha, nothing, decay) + end +end + +function update_penalization!(reg::SDDPRegularization) + reg.ρ *= reg.alpha +end + +function getpenaltyexpr(reg::SDDPRegularization, x, xp) + QuadExpr(reg.ρ*dot(x - xp, x - xp)) +end + +function push_state!(reg::SDDPRegularization, x::Vector, t::Int, k::Int) + incumbents[t+1, k, :] = reg.decay*x + (1-reg.decay)*incumbents[t+1, k, :] +end + +getincumbent(reg::SDDPRegularization, t::Int, k::Int) = reg.incumbents[t+1, k, :] +getavgincumbent(reg::SDDPRegularization, t::Int) = mean(reg.incumbents[t+1, :, :], 2) + diff --git a/src/stoppingtest.jl b/src/stoppingtest.jl index 413da74..0f0688e 100644 --- a/src/stoppingtest.jl +++ b/src/stoppingtest.jl @@ -1,11 +1,9 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# -# Implement the SDDP solver and initializers: -# - functions to initialize value functions -# - functions to build terminal cost +# SDDP stopping criterion ############################################################################# diff --git a/src/utils.jl b/src/utils.jl index 3fffd2b..228accd 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,4 +1,4 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard +# Copyright 2017, V.Leclere, H.Gerard, F.Pacaud, T.Rigaut # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. @@ -7,18 +7,21 @@ # ############################################################################# +import Base: +, show, writecsv """ -Dump Polyhedral functions in a text file. +Write Polyhedral functions in a CSV file. + +$(SIGNATURES) # Arguments -* `dump::String`: +* `filename::String`: Name of output filt * `Vts::Vector{PolyhedralFunction}`: Vector of polyhedral functions to save """ -function dump_polyhedral_functions(dump::AbstractString, Vts::Vector{PolyhedralFunction}) - outfile = open(dump, "w") +function writecsv(filename::AbstractString, Vts::Vector{PolyhedralFunction}) + outfile = open(filename, "w") time = 1 for V in Vts @@ -64,8 +67,7 @@ function read_polyhedral_functions(dump::AbstractString) V[t].betas = vcat(V[t].betas, beta) V[t].numCuts += 1 catch - V[t] = PolyhedralFunction([beta], - reshape(lambda, 1, dim_state), 1) + V[t] = PolyhedralFunction([beta], reshape(lambda, 1, dim_state)) end end return V @@ -90,28 +92,6 @@ function catcuts(Vts::StochDynamicProgramming.PolyhedralFunction...) return StochDynamicProgramming.PolyhedralFunction(betas, lambdas, numcuts) end -""" -Extract a vector stored in a 3D Array - -# Arguments -* `input_array::Array{Float64, 3}`: - array storing the values of vectors -* `nx::Int64`: - Position of vector in first dimension -* `ny::Int64`: - Position of vector in second dimension - -# Return -`Vector{Float64}` -""" -function extract_vector_from_3Dmatrix(input_array::Array{Float64, 3}, - nx::Int64, - ny::Int64) - info("extract_vector_from_3Dmatrix is now deprecated. Use collect instead.") - state_dimension = size(input_array)[3] - return reshape(input_array[nx, ny, :], state_dimension) -end - """ Generate a random state. @@ -129,19 +109,31 @@ end """ -Print in terminal: +Print in stdout: Pass number Upper bound Lower bound exectime + # Arguments +* `io::IO`: * `stats::SDDPStat`: -* `verbose::Int64`: + """ -function print_current_stats(stats::SDDPStat, verbose::Int64) - if (verbose > 0) && (stats.niterations%verbose==0) - print("Pass n\° ", stats.niterations) - (stats.upper_bounds[end] < Inf) && @printf("\tUpper-bound: %.4e", stats.upper_bounds[end]) - @printf("\tLower-bound: %.4e", stats.lower_bounds[end]) - println("\tTime: ", round(stats.exectime[end], 2),"s") - end +function Base.show(io::IO, stats::SDDPStat) + print("Pass n\° ", stats.niterations) + (stats.upper_bounds[end] < Inf) && @printf("\tUpper-bound: %.4e", stats.upper_bounds[end]) + @printf("\tLower-bound: %.4e", stats.lower_bounds[end]) + print("\tTime: ", round(stats.exectime[end], 2),"s") end +"""Check if `k` is congruent with current iteration `it`.""" +checkit(k::Int, it::Int) = k > 0 && it%k == 0 + +function showperformance(stats::SDDPStat) + tbw = sum(stats.solverexectime_bw) + tfw = sum(stats.solverexectime_fw) + titer = sum(stats.exectime) + println("Time in forward pass: $tfw") + println("Time in backward pass: $tbw") + println("Total solver time: $(tfw+tbw)") + println("Total execution time: $titer") +end diff --git a/test/REQUIRE b/test/REQUIRE index b23f66f..df9b20b 100644 --- a/test/REQUIRE +++ b/test/REQUIRE @@ -1,2 +1 @@ -FactCheck Clp diff --git a/test/extensive_formulation.jl b/test/extensive_formulation.jl index ab5d290..c8b7d3a 100644 --- a/test/extensive_formulation.jl +++ b/test/extensive_formulation.jl @@ -1,9 +1,10 @@ ################################################################################ -# Test SDDP functions +# Test extensive formulation ################################################################################ -using FactCheck, StochDynamicProgramming, JuMP, Clp -facts("SDDP algorithm: 1D case") do +using StochDynamicProgramming, Base.Test, Clp + +@testset "Extensive formulation" begin solver = ClpSolver() # SDDP's tolerance: @@ -48,14 +49,14 @@ facts("SDDP algorithm: 1D case") do gap=epsilon, max_iterations=max_iterations) - context("Extensive solving") do + @testset "Extensive solving" begin ef_cost = StochDynamicProgramming.extensive_formulation(model_ef, params)[1] - @fact typeof(ef_cost) --> Float64 + @test isa(ef_cost, Float64) end - context("Unsolvable extensive formulation") do + @testset "Unsolvable extensive formulation" begin x_bounds_ef = [(-2., -1.)] set_state_bounds(model_ef, x_bounds_ef) - @fact_throws extensive_formulation(model_ef, params) + @test_throws ErrorException extensive_formulation(model_ef, params) end end diff --git a/test/prob.jl b/test/prob.jl index b1ca0f9..67c99a0 100644 --- a/test/prob.jl +++ b/test/prob.jl @@ -1,30 +1,30 @@ ################################################################################ # Test probability functions ################################################################################ -using FactCheck, StochDynamicProgramming +using Base.Test, StochDynamicProgramming -facts("Probability functions") do +@testset "Probability functions" begin support = [1, 2, 3] proba = [.2 .5 .3] law = NoiseLaw(support, proba) - @fact typeof(law) --> NoiseLaw - @fact law.supportSize --> 3 + @test typeof(law) == NoiseLaw + @test law.supportSize == 3 dims = (2, 2, 1) scenarios = simulate_scenarios([law, law], 2) - @fact typeof(scenarios) --> Array{Float64, 3} - @fact size(scenarios) --> (2, 2, 1) + @test typeof(scenarios) == Array{Float64, 3} + @test size(scenarios) == (2, 2, 1) # test product of noiselaws: support2 = [4, 5, 6] proba2 = [.3 .3 .4] law2 = NoiseLaw(support2, proba2) law3 = StochDynamicProgramming.noiselaw_product(law, law2) - @fact law3.supportSize --> law.supportSize*law2.supportSize - @fact law3.proba --> vec(proba' * proba2) - @fact size(law3.support)[1] --> size(law.support)[1] + size(law2.support)[1] - @fact law3.support[:, 1] --> [1., 4.] + @test law3.supportSize == law.supportSize*law2.supportSize + @test law3.proba == vec(proba' * proba2) + @test size(law3.support)[1] == size(law.support)[1] + size(law2.support)[1] + @test law3.support[:, 1] == [1., 4.] # Test product of three noiselaws: StochDynamicProgramming.noiselaw_product(law, law2, law) diff --git a/test/runtests.jl b/test/runtests.jl index 3a9490d..c45ae91 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,12 +8,12 @@ using StochDynamicProgramming -using Distributions, Clp, FactCheck, JuMP +using Clp, JuMP +using Base.Test # Test utility functions: include("prob.jl") -include("utils.jl") # Test SDDP: include("sddp.jl") diff --git a/test/sddp.jl b/test/sddp.jl index 2d8b663..27432fe 100644 --- a/test/sddp.jl +++ b/test/sddp.jl @@ -1,10 +1,11 @@ ################################################################################ # Test SDDP functions ################################################################################ -using FactCheck, StochDynamicProgramming, JuMP, Clp +using StochDynamicProgramming, JuMP, Clp +using Base.Test # Test SDDP with a one dimensional stock: -facts("SDDP algorithm: 1D case") do +@testset "SDDP algorithm: 1D case" begin solver = ClpSolver() # SDDP's tolerance: @@ -43,164 +44,114 @@ facts("SDDP algorithm: 1D case") do # Instantiate parameters of SDDP: param = StochDynamicProgramming.SDDPparameters(solver, - passnumber=n_scenarios, - gap=epsilon, - max_iterations=max_iterations, - prune_cuts=0) + passnumber=n_scenarios, + gap=epsilon, + max_iterations=max_iterations, + prune_cuts=0) V = nothing model = StochDynamicProgramming.LinearSPModel(n_stages, u_bounds, - x0, cost, dynamic, laws) + x0, cost, dynamic, laws) set_state_bounds(model, x_bounds) # Test error if bounds are not well specified: - @fact_throws set_state_bounds(model, [(0,1), (0,1)]) + @test_throws ErrorException set_state_bounds(model, [(0,1), (0,1)]) # Generate scenarios for forward simulations: noise_scenarios = simulate_scenarios(model.noises,param.forwardPassNumber) sddp_costs = 0 - context("Linear cost") do + @testset "Linear cost" begin # Compute bellman functions with SDDP: - V, pbs = solve_SDDP(model, param, 0) - @fact typeof(V) --> Vector{StochDynamicProgramming.PolyhedralFunction} - @fact typeof(pbs) --> Vector{JuMP.Model} - @fact length(pbs) --> n_stages - 1 - @fact length(V) --> n_stages - + sddp = solve_SDDP(model, param, 0) + @test isa(sddp, SDDPInterface) + @test typeof(sddp.bellmanfunctions) == Vector{StochDynamicProgramming.PolyhedralFunction} + @test typeof(sddp.solverinterface) == Vector{JuMP.Model} + @test length(sddp.solverinterface) == n_stages - 1 + @test length(sddp.bellmanfunctions) == n_stages + + V = sddp.bellmanfunctions # Test if the first subgradient has the same dimension as state: - @fact size(V[1].lambdas, 2) --> model.dimStates - @fact V[1].numCuts <= n_scenarios*max_iterations + n_scenarios --> true - @fact size(V[1].lambdas, 1) --> V[1].numCuts + @test size(V[1].lambdas, 2) == model.dimStates + @test V[1].numCuts <= n_scenarios*max_iterations + n_scenarios + @test size(V[1].lambdas, 1) == V[1].numCuts # Test upper bounds estimation with Monte-Carlo: n_simulations = 100 - upb = StochDynamicProgramming.estimate_upper_bound(model, param, V, pbs, - n_simulations)[1] - @fact typeof(upb) --> Float64 - - sddp_costs, stocks = forward_simulations(model, param, pbs, noise_scenarios) + upb = StochDynamicProgramming.estimate_upper_bound(model, + param, + V, + sddp.solverinterface, + n_simulations)[1] + @test typeof(upb) == Float64 + + pbs = sddp.solverinterface + sddp_costs, stocks = forward_simulations(model, param, pbs, + noise_scenarios) # Test error if scenarios are not given in the right shape: - @fact_throws forward_simulations(model, param, pbs, [1.]) + @test_throws BoundsError forward_simulations(model, param, pbs, [1.]) # Test computation of optimal control: - aleas = collect(noise_scenarios[1, 1, :]) - opt = StochDynamicProgramming.get_control(model, param, pbs, 1, model.initialState, aleas) - @fact typeof(opt) --> Vector{Float64} + aleas = noise_scenarios[1, 1, :] + opt = StochDynamicProgramming.get_control(model, param, + sddp.solverinterface, + 1, model.initialState, aleas) + @test typeof(opt) == Vector{Float64} # Test display: - StochDynamicProgramming.set_max_iterations(param, 2) param.compute_ub = 0 - V, pbs, stats = solve_SDDP(model, param, V, 2) + sddp = solve_SDDP(model, param, V, 2) end - context("Value functions calculation") do + @testset "Value functions calculation" begin V0 = StochDynamicProgramming.get_lower_bound(model, param, V) + @test isa(V0, Float64) end - context("Hotstart") do + @testset "Hotstart" begin # Test hot start with previously computed value functions: - V, pbs = solve_SDDP(model, param, V, 0) + sddp = solve_SDDP(model, param, V, 0) + @test isa(sddp, SDDPInterface) # Test if costs are roughly the same: - sddp_costs2, stocks = forward_simulations(model, param, pbs, noise_scenarios) - @fact mean(sddp_costs) --> roughly(mean(sddp_costs2)) + sddp_costs2, stocks = forward_simulations(model, param, + sddp.solverinterface, noise_scenarios) + @test mean(sddp_costs) ≈ mean(sddp_costs2) end - context("Cuts pruning") do - v = V[1] - vt = PolyhedralFunction([v.betas[1]; v.betas[1] - 1.], v.lambdas[[1,1],:], 2) - # Check cuts counting: - @fact StochDynamicProgramming.get_total_number_cuts([vt]) --> 2 - - # Check computation of cut value: - @fact StochDynamicProgramming.cutvalue(vt, 1, [0., 0.]) --> v.betas[1] - - # Check computation of optimal cut: - @fact StochDynamicProgramming.optimalcut([0., 0.], vt)[2] --> 1 - - terr = StochDynamicProgramming.ActiveCutsContainer(2) - StochDynamicProgramming.find_level1_cuts!(terr, vt, [0. 0.; 1. 0.]) - @fact terr.numCuts --> 2 - @fact terr.nstates --> 2 - @fact length(terr.territories[1]) --> 2 - @fact length(terr.territories[2]) --> 0 - - # Check heuristic removal: - vt2 = StochDynamicProgramming.level1_cuts_pruning!(model, param, vt, terr) - @fact isa(vt2, StochDynamicProgramming.PolyhedralFunction) --> true - @fact vt2.numCuts --> 1 - @fact vt2.betas[1] --> vt.betas[1] - - # Check exact dominance test: - isactive1 = StochDynamicProgramming.is_cut_relevant(model, 1, vt, param.SOLVER)[1] - isactive2 = StochDynamicProgramming.is_cut_relevant(model, 2, vt, param.SOLVER)[1] - @fact isactive1 --> true - @fact isactive2 --> false - - # Check insertion of pruning algorithms into SDDP solver: - param1 = StochDynamicProgramming.SDDPparameters(solver, - passnumber=n_scenarios, - gap=epsilon, - pruning_algo="exact", - prune_cuts=1, - max_iterations=1) - V1 = solve_SDDP(model, param1, 0)[1] + @testset "Quadratic regularization" begin param2 = StochDynamicProgramming.SDDPparameters(solver, passnumber=n_scenarios, gap=epsilon, - pruning_algo="level1", - prune_cuts=1, - max_iterations=1) - V2 = solve_SDDP(model, param2, 0)[1] - param3 = StochDynamicProgramming.SDDPparameters(solver, - passnumber=n_scenarios, - gap=epsilon, - pruning_algo="exact+", - prune_cuts=1, - max_iterations=1) - V3 = solve_SDDP(model, param3, 0)[1] - - n1 = StochDynamicProgramming.get_total_number_cuts(V1) - n2 = StochDynamicProgramming.get_total_number_cuts(V2) - n3 = StochDynamicProgramming.get_total_number_cuts(V3) - @fact n1 > n2 --> true - @fact n3 > n2 --> true - end - - context("Quadratic regularization") do - param2 = StochDynamicProgramming.SDDPparameters(solver, - passnumber=n_scenarios, - gap=epsilon, - max_iterations=max_iterations, - rho0=1.) + max_iterations=max_iterations) #TODO: fix solver, as Clp cannot solve QP - @fact_throws solve_SDDP(model, param2, 0) + @test_throws ErrorException solve_SDDP(model, param2, 0, + regularization=SDDPRegularization(1., .99)) end # Test definition of final cost with a JuMP.Model: - context("Final cost") do + @testset "Final cost" begin function fcost(model, m) alpha = getvariable(m, :alpha) @constraint(m, alpha == 0.) end # Store final cost in model: model.finalCost = fcost - V, pbs = solve_SDDP(model, param, 0) - V, pbs = solve_SDDP(model, param, V, 0) + solve_SDDP(model, param, 0) + solve_SDDP(model, param, V, 0) end - context("Piecewise linear cost") do + @testset "Piecewise linear cost" begin # Test Piecewise linear costs: model = StochDynamicProgramming.LinearSPModel(n_stages, u_bounds, x0, [cost], dynamic, laws) set_state_bounds(model, x_bounds) - V, pbs = solve_SDDP(model, param, 0) + sddp = solve_SDDP(model, param, 0) end - context("SMIP") do + @testset "SMIP" begin controlCat = [:Bin, :Cont] u_bounds = [(0., 1.), (0., Inf)] model2 = StochDynamicProgramming.LinearSPModel(n_stages, @@ -209,43 +160,45 @@ facts("SDDP algorithm: 1D case") do dynamic, laws, control_cat=controlCat) set_state_bounds(model2, x_bounds) - @fact_throws solve_SDDP(model2, param, 0) + @test_throws ErrorException solve_SDDP(model2, param, 0) end - context("Stopping criterion") do + @testset "Stopping criterion" begin # Compute upper bound every %% iterations: param.compute_ub = 1 - param.maxItNumber = 30 - param.gap = .1 - V, pbs = solve_SDDP(model, param, V, 0) - V0 = StochDynamicProgramming.get_lower_bound(model, param, V) + gap = .1 + sddp = solve_SDDP(model, param, V, 0) + V0 = StochDynamicProgramming.get_lower_bound(model, param, sddp.bellmanfunctions) n_simulations = 1000 - upb = StochDynamicProgramming.estimate_upper_bound(model, param, V, pbs, - n_simulations)[1] - @fact abs((V0 - upb)) < param.gap --> true + upb = StochDynamicProgramming.estimate_upper_bound(model, param, + sddp.bellmanfunctions, + sddp.solverinterface, + n_simulations)[1] + + @test abs((V0 - upb))/V0 < gap end - context("Dump") do + @testset "Dump" begin # Dump V in text file: - StochDynamicProgramming.dump_polyhedral_functions("dump.dat", V) + StochDynamicProgramming.writecsv("dump.dat", V) # Get stored values: Vdump = StochDynamicProgramming.read_polyhedral_functions("dump.dat") - @fact V[1].numCuts --> Vdump[1].numCuts - @fact V[1].betas --> Vdump[1].betas - @fact V[1].lambdas --> Vdump[1].lambdas + @test V[1].numCuts == Vdump[1].numCuts + @test V[1].betas == Vdump[1].betas + @test V[1].lambdas == Vdump[1].lambdas end - context("Compare parameters") do - paramDDP = [param for i in 1:3] - scenarios = StochDynamicProgramming.simulate_scenarios(laws, 1000) - benchmark_parameters(model, paramDDP, scenarios, 12) - end + #= @testset "Compare parameters" begin =# + #= paramDDP = [param for i in 1:3] =# + #= scenarios = StochDynamicProgramming.simulate_scenarios(laws, 1000) =# + #= benchmark_parameters(model, paramDDP, scenarios, 12) =# + #= end =# end # Test SDDP with a two-dimensional stock: -facts("SDDP algorithm: 2D case") do +@testset "SDDP algorithm: 2D case" begin solver = ClpSolver() # SDDP's tolerance: @@ -288,29 +241,33 @@ facts("SDDP algorithm: 2D case") do gap=epsilon, max_iterations=max_iterations) V = nothing - context("Linear cost") do + @testset "Linear cost" begin # Instantiate a SDDP linear model: model = StochDynamicProgramming.LinearSPModel(n_stages, - u_bounds, x0, - cost, - dynamic, laws) + u_bounds, x0, + cost, + dynamic, laws) set_state_bounds(model, x_bounds) # Compute bellman functions with SDDP: - V, pbs, stats = solve_SDDP(model, param, 0) - @fact typeof(V) --> Vector{StochDynamicProgramming.PolyhedralFunction} - @fact typeof(pbs) --> Vector{JuMP.Model} - @fact typeof(stats) --> StochDynamicProgramming.SDDPStat + sddp = solve_SDDP(model, param, 0) + @test isa(sddp, SDDPInterface) + V = sddp.bellmanfunctions + pbs = sddp.solverinterface + stats = sddp.stats + @test typeof(V) == Vector{StochDynamicProgramming.PolyhedralFunction} + @test typeof(pbs) == Vector{JuMP.Model} + @test typeof(stats) == StochDynamicProgramming.SDDPStat # Test if the first subgradient has the same dimension as state: - @fact length(V[1].lambdas[1, :]) --> model.dimStates + @test length(V[1].lambdas[1, :]) == model.dimStates # Test upper bounds estimation with Monte-Carlo: n_simulations = 100 upb = StochDynamicProgramming.estimate_upper_bound(model, param, V, pbs, - n_simulations)[1] - @fact typeof(upb) --> Float64 + n_simulations)[1] + @test typeof(upb) == Float64 # Test a simulation upon given scenarios: @@ -320,22 +277,22 @@ facts("SDDP algorithm: 2D case") do # Compare sddp cost with those given by extensive formulation: ef_cost = StochDynamicProgramming.extensive_formulation(model,param)[1] - @fact typeof(ef_cost) --> Float64 + @test typeof(ef_cost) == Float64 - @fact mean(sddp_costs) --> roughly(ef_cost) + @test mean(sddp_costs) ≈ ef_cost end - context("Dump") do + @testset "Dump" begin # Dump V in text file: - StochDynamicProgramming.dump_polyhedral_functions("dump.dat", V) + StochDynamicProgramming.writecsv("dump.dat", V) # Get stored values: Vdump = StochDynamicProgramming.read_polyhedral_functions("dump.dat") - @fact V[1].numCuts --> Vdump[1].numCuts - @fact V[1].betas --> Vdump[1].betas - @fact V[1].lambdas --> Vdump[1].lambdas + @test V[1].numCuts == Vdump[1].numCuts + @test V[1].betas == Vdump[1].betas + @test V[1].lambdas == Vdump[1].lambdas end end diff --git a/test/sdp.jl b/test/sdp.jl index 54d2e76..257689d 100644 --- a/test/sdp.jl +++ b/test/sdp.jl @@ -1,10 +1,10 @@ ################################################################################ # Test SDDP functions ################################################################################ -using FactCheck, StochDynamicProgramming +using Base.Test, StochDynamicProgramming using StochDynamicProgramming.SDPutils -facts("Indexation for SDP") do +@testset "Indexation for SDP" begin bounds = [(0.1,10.0), (1.2, 4.0), (0.5, 2.0)] steps = [0.1, 0.05, 0.01] @@ -18,17 +18,17 @@ facts("Indexation for SDP") do checkTrue = SDPutils.is_next_state_feasible([0.12,1.3,1.3],3,bounds) - @fact ind --> (4,51,141) - @fact ind2[1] --> roughly(4.2) - @fact ind2[2] --> roughly(52.6) - @fact ind2[3] --> roughly(144.2) - @fact checkFalse --> false - @fact checkTrue --> true + @test ind == (4,51,141) + @test ind2[1] ≈ 4.2 + @test ind2[2] ≈ 52.6 + @test ind2[3] ≈ 144.2 + @test ~checkFalse + @test checkTrue end -facts("SDP algorithm") do +@testset "SDP algorithm" begin # Number of timesteps : TF = 3 @@ -127,7 +127,7 @@ facts("SDP algorithm") do "Exact"); - context("Compare StochDynProgModel constructors") do + @testset "Compare StochDynProgModel constructors" begin modelSDPPiecewise = StochDynamicProgramming.LinearSPModel(TF, @@ -160,20 +160,20 @@ facts("SDP algorithm") do test_costs &= (modelSDPPiecewise.costFunctions[1](t,x,u,w)==convertedSDPmodel.costFunctions(t,x,u,w)) end - @fact test_costs --> true + @test test_costs - @fact convertedSDPmodel.constraints(1,x,u,w) --> true + @test convertedSDPmodel.constraints(1,x,u,w) end - context("Solve and simulate using SDP") do + @testset "Solve and simulate using SDP" begin paramsSDP.infoStructure = "anything" - @fact_throws solve_DP(modelSDP, paramsSDP, false); + @test_throws ErrorException solve_DP(modelSDP, paramsSDP, false); paramsSDP.infoStructure = infoStruct V_sdp = solve_DP(modelSDP, paramsSDP, false); - @fact size(V_sdp) --> (paramsSDP.stateVariablesSizes..., TF) + @test size(V_sdp) == (paramsSDP.stateVariablesSizes..., TF) costs_sdp, stocks_sdp, controls_sdp = StochDynamicProgramming.sdp_forward_single_simulation(modelSDP, paramsSDP, @@ -186,7 +186,7 @@ facts("SDP algorithm") do aleas_scen, V_sdp, true ) - @fact costs_sdp2[1] --> costs_sdp + @test costs_sdp2[1] == costs_sdp x = x0 V_sdp2 = StochDynamicProgramming.sdp_compute_value_functions(modelSDP, paramsSDP, false); @@ -202,8 +202,8 @@ facts("SDP algorithm") do v2 = Vitp2[(1.1,1.1)...] v3 = Vitp3[(1.1,1.1)...] - @fact v1 --> v2 - @fact (v1<=v3) --> true + @test v1 == v2 + @test v1 <= v3 paramsSDP.infoStructure = "DH" costs_sdp3, stocks_sdp3, controls_sdp3 = StochDynamicProgramming.sdp_forward_simulation(modelSDP, @@ -212,7 +212,7 @@ facts("SDP algorithm") do V_sdp3, true ) paramsSDP.infoStructure = "HD" - @fact costs_sdp3[1]>=costs_sdp2[1] --> true + @test costs_sdp3[1] >= costs_sdp2[1] a,b = StochDynamicProgramming.generate_grid(modelSDP, paramsSDP) @@ -222,33 +222,33 @@ facts("SDP algorithm") do u_bounds = modelSDP.ulim u_steps = paramsSDP.controlSteps - @fact length(collect(a)) --> (x_bounds[1][2]-x_bounds[1][1]+x_steps[1])*(x_bounds[2][2]-x_bounds[2][1]+x_steps[2])/(x_steps[1]*x_steps[2]) - @fact length(collect(b)) --> (u_bounds[1][2]-u_bounds[1][1]+u_steps[1])*(u_bounds[2][2]-u_bounds[2][1]+u_steps[2])/(u_steps[1]*u_steps[2]) + @test length(collect(a)) == (x_bounds[1][2]-x_bounds[1][1]+x_steps[1])*(x_bounds[2][2]-x_bounds[2][1]+x_steps[2])/(x_steps[1]*x_steps[2]) + @test length(collect(b)) == (u_bounds[1][2]-u_bounds[1][1]+u_steps[1])*(u_bounds[2][2]-u_bounds[2][1]+u_steps[2])/(u_steps[1]*u_steps[2]) ind = SDPutils.index_from_variable(x, x_bounds, x_steps) - @fact get_bellman_value(modelSDP, paramsSDP, V_sdp2) --> V_sdp2[ind...,1] + @test get_bellman_value(modelSDP, paramsSDP, V_sdp2) == V_sdp2[ind...,1] - @fact size(V_sdp) --> (paramsSDP.stateVariablesSizes..., TF) - @fact V_sdp2[1,1,1] <= V_sdp3[1,1,1] --> true + @test size(V_sdp) == (paramsSDP.stateVariablesSizes..., TF) + @test V_sdp2[1,1,1] <= V_sdp3[1,1,1] - @fact size(stocks_sdp) --> (3,1,2) - @fact size(controls_sdp) --> (2,1,2) + @test size(stocks_sdp) == (3,1,2) + @test size(controls_sdp) == (2,1,2) state_ref = zeros(2) state_ref[1] = stocks_sdp[2,1,1] state_ref[2] = stocks_sdp[2,1,2] w = [4] - @fact_throws get_control(modelSDP,paramsSDP,V_sdp3, 1, x) - @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x, w)[1] >= CONTROL_MIN) --> true - @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x, w)[1] <= CONTROL_MAX) --> true + @test_throws ErrorException get_control(modelSDP,paramsSDP,V_sdp3, 1, x) + @test (get_control(modelSDP,paramsSDP,V_sdp3, 1, x, w)[1] >= CONTROL_MIN) + @test (get_control(modelSDP,paramsSDP,V_sdp3, 1, x, w)[1] <= CONTROL_MAX) paramsSDP.infoStructure = "DH" - @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x)[1] >= CONTROL_MIN) --> true - @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x)[1] <= CONTROL_MAX) --> true + @test (get_control(modelSDP,paramsSDP,V_sdp3, 1, x)[1] >= CONTROL_MIN) + @test (get_control(modelSDP,paramsSDP,V_sdp3, 1, x)[1] <= CONTROL_MAX) - @fact size(stocks_sdp) --> (3,1,2) - @fact size(controls_sdp) --> (2,1,2) + @test size(stocks_sdp) == (3,1,2) + @test size(controls_sdp) == (2,1,2) end