diff --git a/.travis.yml b/.travis.yml index 5321b7a..1afb4b2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -3,6 +3,7 @@ os: - linux julia: - 0.5 + - 0.6 notifications: email: false diff --git a/README.md b/README.md index 43d9bb4..71673b4 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ **WARNING:** *This package is currently in development. Any help or feedback is appreciated.* -**Latest release:** 0.3.0 +**Latest release:** 0.4.0 | **Documentation** | **Build Status** | **Social** | |:-----------------:|:----------------:|:----------:| @@ -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,7 +54,10 @@ control strategies. ## Installation -Installing StochDynamicProgramming is an easy process. Open Julia and enter +Installing StochDynamicProgramming is an easy process. +Currently, the package depends upon `StochasticDualDynamicProgramming.jl`, which is not +yet registered in Julia's METADATA. To install the package, +open Julia and enter ```julia julia> Pkg.update() diff --git a/REQUIRE b/REQUIRE index f20beb6..eb96af6 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,6 +1,7 @@ julia 0.5 -JuMP 0.16 +JuMP 0.17 Distributions ProgressMeter Interpolations -Iterators +CutPruners 0.0.2 +Compat 0.18 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/doc/sddp_api.rst b/doc/sddp_api.rst index 0b02161..ba606ed 100644 --- a/doc/sddp_api.rst +++ b/doc/sddp_api.rst @@ -20,8 +20,8 @@ but allows to use either piecewise linear or quadratic costs. To define a `LinearSPModel`, the constructor is:: spmodel = LinearSPModel( - nstage, # number of stages - ubounds, # bounds of control + n_stage, # number of stages + u_bounds, # bounds of control x0, # initial state cost, # cost function dynamic, # dynamic @@ -36,8 +36,8 @@ Default parameters ^^^^^^^^^^^^^^^^^^ You should at least specify these parameters to define a `LinearSPModel`: -- `nstage` (Int): number of stages in the stochastic multistage problem -- `ubounds` (list of tuple): bounds upon control, defined as a sequence of tuple :code:`(umin, umax)`. +- `n_stage` (Int): number of stages in the stochastic multistage problem +- `u_bounds` (list of tuple): bounds upon control, defined as a sequence of tuple :code:`(umin, umax)`. - `x0` (`Vec{Float64}`): initial state - `cost` (`Function`): cost function as a function of time, state, control and noise returning a Float - `dynamic` (`Function`): system's dynamic as a function of time, state, control and noise returning a vector diff --git a/examples/battery_storage_parallel.jl b/examples/battery_storage_parallel.jl index 1bbed4c..40e471f 100644 --- a/examples/battery_storage_parallel.jl +++ b/examples/battery_storage_parallel.jl @@ -106,10 +106,10 @@ println("library loaded") controlSteps, infoStruct) end -Vs = StochDynamicProgramming.solve_DP(spmodel,paramSDP, 1) +Vs = StochDynamicProgramming.solve_dp(spmodel,paramSDP, 1) lb_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs) println("Value obtained by SDP: "*string(lb_sdp)) -costsdp, states, stocks = StochDynamicProgramming.sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs) +costsdp, states, stocks = StochDynamicProgramming.forward_simulations(spmodel,paramSDP,Vs,scenarios) println(mean(costsdp)) 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..8bbea04 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 @@ -67,7 +66,7 @@ function solve_determinist_problem() @variable(m, 0. <= u[1:N_STAGES-1] <= 7) @variable(m, 0. <= s[1:N_STAGES-1] <= 7) - @objective(m, Min, sum{COST[i]*u[i], i = 1:N_STAGES-1}) + @objective(m, Min, sum(COST[i]*u[i] for i = 1:N_STAGES-1)) for i in 1:(N_STAGES-1) @constraint(m, x[i+1] - x[i] + u[i] + s[i] - alea_year[i] == 0) @@ -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..5e51e1c 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") ################################################## @@ -58,11 +60,11 @@ end function final_cost_dams(model, m) # Here, model is the optimization problem at time T - 1 # so that xf (x future) is the final stock - alpha = JuMP.getvariable(m, :alpha) - w = JuMP.getvariable(m, :w) - x = JuMP.getvariable(m, :x) - u = JuMP.getvariable(m, :u) - xf = JuMP.getvariable(m, :xf) + alpha = m[:alpha] + w = m[:w] + x = m[:x] + u = m[:u] + xf = m[:xf] @JuMP.variable(m, z1 >= 0) @JuMP.variable(m, z2 >= 0) @JuMP.variable(m, z3 >= 0) @@ -115,10 +117,8 @@ function init_problem() # Add bounds for stocks: 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) - params = SDDPparameters(solver, + params = SDDPparameters(SOLVER, passnumber=FORWARD_PASS, compute_ub=10, gap=EPSILON, @@ -128,5 +128,5 @@ end # Solve the problem: model, params = init_problem() -V, pbs = @time solve_SDDP(model, params, 1) +sddp = @time solve_SDDP(model, params, 2, 1) diff --git a/examples/multistock-example.jl b/examples/multistock-example.jl index a08b6e5..a08ec98 100644 --- a/examples/multistock-example.jl +++ b/examples/multistock-example.jl @@ -16,8 +16,8 @@ using StochDynamicProgramming, Clp println("library loaded") run_sddp = true # false if you don't want to run sddp -run_sdp = true # false if you don't want to run sdp -test_simulation = true # false if you don't want to test your strategies +run_sdp = false # false if you don't want to run sdp +test_simulation = false # false if you don't want to test your strategies ######## Optimization parameters ######## # choose the LP solver used. @@ -25,7 +25,7 @@ const SOLVER = ClpSolver() # require "using Clp" #const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0) # require "using CPLEX" # convergence test -const MAX_ITER = 10 # number of iterations of SDDP +const MAX_ITER = 100 # number of iterations of SDDP const step = 0.1 # discretization step of SDP @@ -81,10 +81,11 @@ if run_sddp println("Starting resolution by SDDP") # 10 forward pass, stop at MAX_ITER paramSDDP = SDDPparameters(SOLVER, - passnumber=10, + passnumber=1, 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 = @time solve_SDDP(spmodel, paramSDDP, 2, # display information every 2 iterations + stopcrit=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 @@ -101,7 +102,7 @@ if run_sdp spmodel_sdp = StochDynamicProgramming.build_sdpmodel_from_spmodel(spmodel) spmodel_sdp.constraints = constraints_dp - Vs = solve_DP(spmodel_sdp, paramSDP, 1) + Vs = solve_dp(spmodel_sdp, paramSDP, 1) value_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs) println("Value obtained by SDP: "*string(round(value_sdp,4))) toc(); println(); @@ -111,8 +112,8 @@ 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) - costsdp, states, controls = sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs) + costsddp, stocks = forward_simulations(spmodel, paramSDDP, sddp.solverinterface, scenarios) + costsdp, states, controls = forward_simulations(spmodel,paramSDP, Vs, scenarios) println("Simulated relative gain of sddp over sdp: " *string(round(200*mean(costsdp-costsddp)/abs(mean(costsddp+costsdp)),3))*"%") end diff --git a/examples/stock-example.jl b/examples/stock-example.jl index 1f9dd88..061645d 100644 --- a/examples/stock-example.jl +++ b/examples/stock-example.jl @@ -21,26 +21,26 @@ test_simulation = true # false if you don't want to test your strategies ######## Optimization parameters ######## # choose the LP solver used. -const SOLVER = ClpSolver() # require "using Clp" +SOLVER = ClpSolver() # require "using Clp" #const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0) # require "using CPLEX" # convergence test -const MAX_ITER = 10 # number of iterations of SDDP -const step = 0.1 # discretization step of SDP +MAX_ITER = 10 # number of iterations of SDDP +step = 0.01 # discretization step of SDP ######## Stochastic Model Parameters ######## -const N_STAGES = 6 # number of stages of the SP problem -const COSTS = [sin(3*t)-1 for t in 1:N_STAGES] +N_STAGES = 6 # number of stages of the SP problem +COSTS = [sin(3*t)-1 for t in 1:N_STAGES-1] #const COSTS = rand(N_STAGES) # randomly generating deterministic costs -const CONTROL_MAX = 0.5 # bounds on the control -const CONTROL_MIN = 0 +CONTROL_MAX = 0.5 # bounds on the control +CONTROL_MIN = 0 -const XI_MAX = 0.3 # bounds on the noise -const XI_MIN = 0 -const N_XI = 10 # discretization of the noise +XI_MAX = 0.3 # bounds on the noise +XI_MIN = 0 +N_XI = 10 # discretization of the noise -const S0 = 0.5 # initial stock +S0 = 0.5 # initial stock # create law of noises proba = 1/N_XI*ones(N_XI) # uniform probabilities @@ -73,8 +73,8 @@ 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 + lb_sddp = StochDynamicProgramming.get_lower_bound(spmodel, paramSDDP, sddp.bellmanfunctions) println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4))) toc(); println(); end @@ -87,7 +87,7 @@ if run_sdp controlSteps = [step] # discretization step of the control infoStruct = "HD" # noise at time t is known before taking the decision at time t paramSDP = SDPparameters(spmodel, stateSteps, controlSteps, infoStruct) - Vs = solve_DP(spmodel,paramSDP, 1) + Vs = solve_dp(spmodel,paramSDP, 1) value_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs) println("Value obtained by SDP: "*string(round(value_sdp,4))) toc(); println(); @@ -109,8 +109,8 @@ 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) - costsdp, states, controls = sdp_forward_simulation(spmodel,paramSDP,scenarios,Vs) + costsddp, stocks = forward_simulations(spmodel, paramSDDP, sddp.solverinterface, scenarios) + costsdp, states, controls = forward_simulations(spmodel,paramSDP, Vs, scenarios) println("Simulated relative gain of sddp over sdp: " *string(round(200*mean(costsdp-costsddp)/abs(mean(costsddp+costsdp)),3))*"%") end diff --git a/src/SDDPoptimize.jl b/src/SDDPoptimize.jl index 65928d1..d8fad1f 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 @@ -21,34 +25,38 @@ fulfilled. the stochastic problem we want to optimize * `param::SDDPparameters`: the parameters of the SDDP algorithm -* `verbose::Int64`: - Default is `0` - If non null, display progression in terminal every - `n` iterations, where `n` is the number specified by display. +* `verbosity::Int64`: + Default is `0`, higher gives more printed information +* `verbose_it::Int64`: + Default is `1` + If verbosity >1 , display progression in terminal every + `verbose_it` iterations. # 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, verbosity=0::Int64, verbose_it=1::Int64; + stopcrit::AbstractStoppingCriterion=IterLimit(param.max_iterations), + prunalgo::AbstractCutPruningAlgo=CutPruners.AvgCutPruningAlgo(-1), + regularization=nothing) + sddp = SDDPInterface(model, param, + stopcrit, + prunalgo, + verbosity=verbosity, + verbose_it=verbose_it, + regularization=regularization) # Run SDDP: - sddp_stats = run_SDDP!(model, param, V, problems, verbose) - return V, problems, sddp_stats + 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. @@ -60,127 +68,132 @@ fulfilled. the parameters of the SDDP algorithm * `V::Vector{PolyhedralFunction}`: current lower approximation of Bellman functions -* `verbose::Int64`: - Default is `0` - If non null, display progression in terminal every - `n` iterations, where `n` is the number specified by display. +* `verbosity::Int64`: + Default is `0`, higher gives more printed information +* `verbose_it::Int64`: + Default is `1` + If verbosity >1 , display progression in terminal every + `verbose_it` iterations. # 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}, verbosity=0::Int64, verbose_it=1::Int64; + stopcrit::AbstractStoppingCriterion=IterLimit(param.max_iterations), + prunalgo::AbstractCutPruningAlgo=CutPruners.AvgCutPruningAlgo(-1)) + + sddp = SDDPInterface(model, param, + stopcrit, + prunalgo, V, + verbosity=verbosity, + verbose_it=verbose_it) + 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) + (sddp.verbosity > 0) && println("Initialize cuts") + 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) + (sddp.verbosity > 0) && println("Starting SDDP iterations") + while !stop(sddp.stopcrit, stats, stats) # Time execution of current pass: tic() #################### # Forward pass : compute stockTrajectories - costs, stockTrajectories, callsolver_forward = forward_pass!(model, param, V, problems) + (sddp.verbosity > 2) && checkit(sddp.verbose_it, sddp.stats.niterations) && println("start forward pass") + costs, trajectories = forward_pass!(sddp) #################### # Backward pass : update polyhedral approximation of Bellman functions - callsolver_backward = backward_pass!(model, param, V, problems, stockTrajectories, model.noises) + (sddp.verbosity > 2) && checkit(sddp.verbose_it, sddp.stats.niterations) && println("start backward pass") + backward_pass!(sddp, trajectories) #################### # 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 + (param.prune) && prune!(sddp, trajectories) + + #################### + # In iteration lower bound estimation + lwb = lowerbound(sddp) #################### # In iteration upper bound estimation - upb = in_iteration_upb_estimation(model, param, stats.niterations+1, verbose, - upperbound_scenarios, upb, problems) + (sddp.verbosity > 2) && checkit(sddp.verbose_it, sddp.stats.niterations) && println("start in-iteration upperbound estimation") + upb = in_iteration_upb_estimation(model, param, stats.niterations+1, sddp.verbosity, + upperbound_scenarios, upb, + sddp.solverinterface) - updateSDDPStat!(stats, callsolver_forward+callsolver_backward, lwb, upb, time_pass) - print_current_stats(stats,verbose) + updateSDDP!(sddp, lwb, upb, time_pass, trajectories) - #################### - # Stopping test - stopping_test = test_stopping_criterion(param,stats) + (sddp.verbosity > 1) && checkit(sddp.verbose_it, 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 end """Display final results once SDDP iterations are finished.""" -function display_final_solution(model::SPModel, param::SDDPparameters, V, - problems, stats::SDDPStat, verbose::Int64) - if (verbose>0) && (param.compute_ub >= 0) - lwb = get_bellman_value(model, param, 1, V[1], model.initialState) +function finalpass!(sddp::SDDPInterface) + model = sddp.spmodel + param = sddp.params + V = sddp.bellmanfunctions + problems = sddp.solverinterface + stats = sddp.stats + + if (sddp.verbosity>0) && (param.compute_ub >= 0) + lwb = lowerbound(sddp) if (param.compute_ub == 0) || (param.monteCarloSize > 0) - (verbose > 0) && println("Compute final upper-bound with ", + println("Compute final upper-bound with ", param.monteCarloSize, " scenarios...") upb, σ, tol = estimate_upper_bound(model, param, V, problems, param.monteCarloSize) else @@ -189,7 +202,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,40 +210,81 @@ 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.params.prune) && sync!(sddp) + (sddp.verbosity >2 )&& println("Reloading JuMP model") + sddp.solverinterface = hotstart_SDDP(sddp.spmodel, + sddp.params, + sddp.bellmanfunctions) + end + + # Update regularization + if !isnull(sddp.regularizer) + (sddp.verbosity >3) && println("Updating regularization ") + 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 +* `verbosity::Int64`: + Default is `0`, higher gives more printed information """ -function build_terminal_cost!(model::SPModel, problem::JuMP.Model, Vt::PolyhedralFunction) +function build_terminal_cost!(model::SPModel, + problem::JuMP.Model, + Vt::PolyhedralFunction, + verbosity::Int64=0) # if shape is PolyhedralFunction, build terminal cost with it: - alpha = getvariable(problem, :alpha) - xf = getvariable(problem, :xf) + alpha = problem[:alpha] + xf = problem[:xf] t = model.stageNumber -1 if isa(Vt, PolyhedralFunction) + (verbosity >3) && println("Building final cost") for i in 1:Vt.numCuts lambda = vec(Vt.lambdas[i, :]) - @constraint(problem, Vt.betas[i] + dot(lambda, xf) <= alpha) + if model.info == :HD + @constraint(problem, Vt.betas[i] + dot(lambda, xf) <= alpha) + elseif model.info == :DH + for ww=1:length(model.noises[t].proba) + @constraint(problem, Vt.betas[i] + dot(lambda, xf[:, ww]) <= alpha[ww]) + end + end end else - @constraint(problem, alpha >= 0) + # else, by default terminal cost is equal to 0 + @constraint(problem, alpha .>= 0) end 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 @@ -246,10 +300,15 @@ linear problem. * `Array::JuMP.Model`: """ function build_models(model::SPModel, param::SDDPparameters) - return JuMP.Model[build_model(model, param, t) for t=1:model.stageNumber-1] + if model.info == :HD + return JuMP.Model[build_model(model, param, t) for t=1:model.stageNumber-1] + else + return JuMP.Model[build_model_dh(model, param, t) for t=1:model.stageNumber-1] + end end -function build_model(model, param, t) + +function build_model(model, param, t,verbosity::Int64=0) m = Model(solver=param.SOLVER) nx = model.dimStates @@ -268,11 +327,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) @@ -293,18 +352,67 @@ function build_model(model, param, t) @objective(m, Min, cost + alpha) end + # store number of cuts + m.ext[:ncuts] = 0 + # Add binary variable if problem is a SMIP: if model.IS_SMIP m.colCat[2*nx+1:2*nx+nu] = model.controlCat end + (verbosity >5) && print(m) return m end +"""Build model in Decision-Hazard.""" +function build_model_dh(model, param, t, verbosity::Int64=0) + m = Model(solver=param.SOLVER) + law = model.noises + + nx = model.dimStates + nu = model.dimControls + nw = model.dimNoises + + ns = law[t].supportSize + ξ = collect(law[t].support[:, :]) + πp = law[t].proba + + @variable(m, model.xlim[i][1] <= x[i=1:nx] <= model.xlim[i][2]) + @variable(m, model.ulim[i][1] <= u[i=1:nu] <= model.ulim[i][2]) + @variable(m, model.xlim[i][1] <= xf[i=1:nx, j=1:ns]<= model.xlim[i][2]) + @variable(m, alpha[1:ns]) + + m.ext[:cons] = @constraint(m, state_constraint, x .== 0) + + for j=1:ns + @constraint(m, xf[:, j] .== model.dynamics(t, x, u, ξ[:, j])) + end + + # add objective as minimization of expectancy: + try + @objective(m, Min, + sum(πp[j]*(model.costFunctions(t, x, u, ξ[:, j]) + + alpha[j]) for j in 1:ns)) + catch + @objective(m, Min, + sum(πp[j]*(model.costFunctions(m, t, x, u, ξ[:, j]) + + alpha[j]) for j in 1:ns)) + end + + # store number of cuts + m.ext[:ncuts] = 0 + + (verbosity >5) && print(m) + return m +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 +432,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 +443,34 @@ 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) end @@ -356,6 +478,8 @@ end Initialize JuMP.Model vector with a previously computed PolyhedralFunction vector. +$(SIGNATURES) + # Arguments * `model::SPModel`: Parametrization of the problem @@ -387,7 +511,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 +546,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 +584,8 @@ end """ Compute optimal control at point xt and time t. +$(SIGNATURES) + # Arguments * `model::SPModel`: Parametrization of the problem @@ -463,13 +605,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 @@ -481,12 +625,19 @@ Add several cuts to JuMP.Model from a PolyhedralFunction Cuts are stored in V """ function add_cuts_to_model!(model::SPModel, t::Int64, problem::JuMP.Model, V::PolyhedralFunction) - alpha = getvariable(problem, :alpha) - xf = getvariable(problem, :xf) + alpha = problem[:alpha] + xf = problem[:xf] for i in 1:V.numCuts lambda = vec(V.lambdas[i, :]) - @constraint(problem, V.betas[i] + dot(lambda, xf) <= alpha) + if model.info == :HD + @constraint(problem, V.betas[i] + dot(lambda, xf) <= alpha) + elseif model.info == :DH + for j in 1:model.noises[t].supportSize + @constraint(problem, V.betas[i] + dot(lambda, xf[:, j]) <= alpha[j]) + end + end end + problem.ext[:ncuts] = V.numCuts end diff --git a/src/SDPoptimize.jl b/src/SDPoptimize.jl deleted file mode 100644 index 584c107..0000000 --- a/src/SDPoptimize.jl +++ /dev/null @@ -1,667 +0,0 @@ -# Copyright 2015, Vincent Leclere, Francois Pacaud, Henri Gerard and -# Tristan 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 -using Iterators -using Interpolations - - -""" -Compute interpolation of the value function at time t - -# Arguments -* `model::SPmodel`: -* `dim_states::Int`: - the number of state variables -* `v::Array`: - the value function to interpolate -* `time::Int`: - time at which we have to interpolate V - -# Return -* Interpolation - the interpolated value function (working as an array with float indexes) - -""" -function value_function_interpolation( dim_states, V, time::Int) - return interpolate(V[[Colon() for i in 1:dim_states]...,time], BSpline(Linear()), OnGrid()) -end - - -""" -Compute the cartesian products of discretized state and control spaces - -# Arguments -* `model::SPmodel`: - the model of the problem -* `param::SDPparameters`: - the parameters of the problem - -# Return -* Iterators: product_states and product_controls - the cartesian product iterators for both states and controls - -""" -function generate_grid(model::SPModel, param::SDPparameters) - product_states = product([model.xlim[i][1]:param.stateSteps[i]:model.xlim[i][2] for i in 1:model.dimStates]...) - product_controls = product([model.ulim[i][1]:param.controlSteps[i]:model.ulim[i][2] for i in 1:model.dimControls]...) - return product_states, product_controls -end - - -""" -Transform a general SPmodel into a StochDynProgModel - -# Arguments -* `model::SPmodel`: - the model of the problem -* `param::SDPparameters`: - the parameters of the problem - -# Return -* `sdpmodel::StochDynProgModel: - the corresponding StochDynProgModel - -""" -function build_sdpmodel_from_spmodel(model::SPModel) - - function zero_fun(x) - return 0 - end - - if isa(model,LinearSPModel) - function cons_fun(t,x,u,w) - return true - end - if in(:finalCostFunction,fieldnames(model)) - SDPmodel = StochDynProgModel(model, model.finalCostFunction, cons_fun) - else - SDPmodel = StochDynProgModel(model, zero_fun, cons_fun) - end - elseif isa(model,StochDynProgModel) - SDPmodel = model - else - error("cannot build StochDynProgModel from current SPmodel. You need to implement - a new StochDynProgModel constructor.") - end - - return SDPmodel -end - - -""" -Dynamic programming algorithm to compute optimal value functions -by backward induction using bellman equation in the finite horizon case. -The information structure can be Decision Hazard (DH) or Hazard Decision (HD) - -# Arguments -* `model::SPmodel`: - the DPSPmodel of our problem -* `param::SDPparameters`: - the parameters for the SDP algorithm -* `display::Int`: - the output display or verbosity parameter - -# Return -* `value_functions::Array`: - the vector representing the value functions as functions of the state - of the system at each time step - -""" -function solve_DP(model::SPModel, param::SDPparameters, display=0::Int64) - SDPmodel = build_sdpmodel_from_spmodel(model::SPModel) - # Start of the algorithm - V = sdp_compute_value_functions(SDPmodel, param, display) - return V -end - - -""" -Compute the value function at time t using bellman equation -and knowing value function at time t+1 - -# Parameters -* `sampling::iz`: (Int) - number of noises samples (number of outcomes if the probability laws are discrete, - number of monte carlo samples otherwise) -* `samples::Array{Float64}`: - arrays of the noises samples/realizations -* `probas::Array{Float64}`: - array of probabilities of samples -* `u_bounds::Array{Tuple{Float64}}`: - array of lower and upper bounds of controls -* `x_bounds::Array{Tuple{Float64}}`: - array of lower and upper bounds of states -* `x_steps::Array{Float64}`: - array of discretization steps for states space -* `x_dim::Int`: - number of state variables -* `product_states::Array{Float64}`: - discretized state space -* `product_controls::Array{Float64}`: - discretized control space -* `dynamics::Function`: - dynamics function of the time step, state, control and randomness returning next state -* `contraints::Function`: - constraints function of the time step, state, control and randomness returning boolean -* `cost::Function`: - cost function of the time step, state, control and randomness returning the instantaneous cost -* `V::Array or SharedArray`: - the array containing the discretized value functions at each time step -* `Vitp::Interpolations`: - the interpolated value functions -* `t::Float64`: - the time step -* `info_struc::String`: - the information structure "HD" for hazard-decision - or "DH" for decision-hazard - -""" -function compute_V_given_t(sampling_size, samples, probas, u_bounds, x_bounds, - x_steps, x_dim, product_states, product_controls, - dynamics, constraints, cost, V, Vitp, t, info_struc) - - if info_struc == "DH" - @sync @parallel for indx in 1:length(product_states) - SDPutils.compute_V_given_x_t_DH(sampling_size, samples, - probas, u_bounds, x_bounds, - x_steps, x_dim, product_controls, - dynamics, constraints, cost, V, Vitp, - t, product_states[indx]) - end - elseif info_struc == "HD" - @sync @parallel for indx in 1:length(product_states) - SDPutils.compute_V_given_x_t_HD(sampling_size, samples, probas, - u_bounds, x_bounds, x_steps, x_dim, - product_controls, dynamics, - constraints, cost, V, Vitp, - t, product_states[indx]) - end - else - error("Information structure should be HD or DH") - end -end - - -""" -Dynamic Programming algorithm to compute optimal value functions - -# Parameters -* `model::StochDynProgModel`: - the StochDynProgModel of the problem -* `param::SDPparameters`: - the parameters for the algorithm -* `display::Int`: - the output display or verbosity parameter - -# Returns -* `value_functions::Array`: - the vector representing the value functions as functions of the state - of the system at each time step - -""" -function sdp_compute_value_functions(model::StochDynProgModel, - param::SDPparameters, - display=0::Int64) - - TF = model.stageNumber - next_state = zeros(Float64, model.dimStates) - - u_bounds = model.ulim - x_bounds = model.xlim - x_steps = param.stateSteps - x_dim = model.dimStates - - dynamics = model.dynamics - constraints = model.constraints - cost = model.costFunctions - - law = model.noises - - #Compute cartesian product spaces - product_states, product_controls = generate_grid(model, param) - - product_states = collect(product_states) - product_controls = collect(product_controls) - - V = SharedArray{Float64}(zeros(Float64, param.stateVariablesSizes..., TF)) - - #Compute final value functions - for x in product_states - ind_x = SDPutils.index_from_variable(x, x_bounds, x_steps) - V[ind_x..., TF] = model.finalCostFunction(x) - end - - if param.expectation_computation!="MonteCarlo" && param.expectation_computation!="Exact" - warn("param.expectation_computation should be 'MonteCarlo' or 'Exact'. Defaulted to 'exact'") - param.expectation_computation="Exact" - end - - #Construct a progress meter - p = 0 - if display > 0 - p = Progress((TF-1), 1) - println("[SDP] Starting value functions computation:") - end - - # Loop over time: - for t = (TF-1):-1:1 - - if display > 0 - next!(p) - end - - if (param.expectation_computation=="MonteCarlo") - sampling_size = param.monteCarloSize - samples = [sampling(law,t) for i in 1:sampling_size] - probas = (1/sampling_size) - else - sampling_size = law[t].supportSize - samples = law[t].support - probas = law[t].proba - end - - Vitp = value_function_interpolation(x_dim, V, t+1) - - compute_V_given_t(sampling_size, samples, probas, u_bounds, x_bounds, - x_steps, x_dim, product_states, product_controls, - dynamics, constraints, cost, V, Vitp, t - , param.infoStructure) - - end - return V -end - - -""" -Get the optimal value of the problem from the optimal Bellman Function - -# Arguments -* `model::SPmodel`: - the DPSPmodel of our problem -* `param::SDPparameters`: - the parameters for the SDP algorithm -* `V::Array{Float64}`: - the Bellman Functions - -# Return -* `V_x0::Float64`: - -""" -function get_bellman_value(model::SPModel, param::SDPparameters, V) - ind_x0 = SDPutils.real_index_from_variable(model.initialState, model.xlim, param.stateSteps) - Vi = value_function_interpolation(model.dimStates, V, 1) - return Vi[ind_x0...,1] -end - - -""" -Simulation of optimal trajectories given model and Bellman functions - -# Arguments -* `model::SPmodel`: - the SPmodel of our problem -* `param::SDPparameters`: - the parameters for the SDP algorithm -* `scenarios::Array`: - the scenarios of uncertainties realizations we want to simulate on - scenarios[t,k,:] is the alea at time t for scenario k -* `V::Array`: - the vector representing the value functions as functions of the state - of the system at each time step -* `display::Bool`: - the output display or verbosity parameter - -# Return -* `costs::Vector{Float64}`: - the cost of the optimal control over the scenario provided -* `states::Array{Float64}`: - the state of the controlled system at each time step -* `controls::Array{Float64}`: - the controls applied to the system at each time step - -""" -function sdp_forward_simulation(model::SPModel, param::SDPparameters, - scenarios::Array{Float64,3}, - V, - display=false::Bool) - - SDPmodel = build_sdpmodel_from_spmodel(model) - TF = SDPmodel.stageNumber - nb_scenarios = size(scenarios)[2] - - costs = zeros(nb_scenarios) - states = zeros(TF,nb_scenarios,model.dimStates) - controls = zeros(TF-1,nb_scenarios,model.dimControls) - - - for k = 1:nb_scenarios - costs[k], states[:,k,:], controls[:,k,:] = sdp_forward_single_simulation(SDPmodel, - param,scenarios[:,k,:],model.initialState,V,display) - end - - return costs, states, controls -end - - -""" -Get the optimal control at time t knowing the state of the system in the decision hazard case - -# Arguments -* `model::SPmodel`: - the DPSPmodel of our problem -* `param::SDPparameters`: - the parameters for the SDP algorithm -* `V::Array{Float64}`: - the Bellman Functions -* `t::int`: - the time step -* `x::Array`: - the state variable -* `w::Array`: - the alea realization - -# Return -* `V_x0::Float64`: - -""" -function get_control(model::SPModel,param::SDPparameters,V, t::Int64, x::Array) - - if(param.infoStructure != "DH") - error("Infostructure must be decision-hazard.") - end - SDPmodel = build_sdpmodel_from_spmodel(model) - - product_controls = product([SDPmodel.ulim[i][1]:param.controlSteps[i]:SDPmodel.ulim[i][2] for i in 1:SDPmodel.dimControls]...) - - law = SDPmodel.noises - best_control = tuple() - Vitp = value_function_interpolation(SDPmodel.dimStates, V, t+1) - - u_bounds = SDPmodel.ulim - x_bounds = SDPmodel.xlim - x_steps = param.stateSteps - - best_V = Inf - - for u in product_controls - - count_admissible_w = 0. - current_V = 0. - - if (param.expectation_computation=="MonteCarlo") - sampling_size = param.monteCarloSize - samples = [sampling(law,t) for i in 1:sampling_size] - probas = (1/sampling_size) - else - sampling_size = law[t].supportSize - samples = law[t].support - probas = law[t].proba - end - - for w = 1:sampling_size - - w_sample = samples[:, w] - proba = probas[w] - - next_state = SDPmodel.dynamics(t, x, u, w_sample) - - if SDPmodel.constraints(t, x, u, w_sample)&&SDPutils.is_next_state_feasible(next_state, model.dimStates, model.xlim) - ind_next_state = SDPutils.real_index_from_variable(next_state, x_bounds, x_steps) - next_V = Vitp[ind_next_state...] - current_V += proba *(SDPmodel.costFunctions(t, x, u, w_sample) + next_V) - count_admissible_w = count_admissible_w + proba - end - end - current_V = current_V/count_admissible_w - if (current_V < best_V)&(count_admissible_w>0) - best_control = u - best_V = current_V - end - end - - return best_control -end - - -""" -Get the optimal control at time t knowing the state of the system and the alea in the hazard decision case - -# Arguments -* `model::SPmodel`: - the DPSPmodel of our problem -* `param::SDPparameters`: - the parameters for the SDP algorithm -* `V::Array{Float64}`: - the Bellman Functions -* `t::int`: - the time step -* `x::Array`: - the state variable -* `w::Array`: - the alea realization - -# Return -* optimal control (tuple(Float64)) - -""" -function get_control(model::SPModel,param::SDPparameters,V, t::Int64, x::Array, w::Array) - - if(param.infoStructure != "HD") - error("Infostructure must be hazard-decision.") - end - - SDPmodel = build_sdpmodel_from_spmodel(model) - - product_controls = product([SDPmodel.ulim[i][1]:param.controlSteps[i]:SDPmodel.ulim[i][2] for i in 1:SDPmodel.dimControls]...) - - law = SDPmodel.noises - best_control = tuple() - Vitp = value_function_interpolation(SDPmodel.dimStates, V, t+1) - - u_bounds = SDPmodel.ulim - x_bounds = SDPmodel.xlim - x_steps = param.stateSteps - - best_V = Inf - - for u = product_controls - - next_state = SDPmodel.dynamics(t, x, u, w) - - if SDPmodel.constraints(t, x, u, w)&&SDPutils.is_next_state_feasible(next_state, model.dimStates, model.xlim) - ind_next_state = SDPutils.real_index_from_variable(next_state, x_bounds, x_steps) - next_V = Vitp[ind_next_state...] - current_V = SDPmodel.costFunctions(t, x, u, w) + next_V - if (current_V < best_V) - best_control = u - best_state = SDPmodel.dynamics(t, x, u, w) - best_V = current_V - end - end - - end - - return best_control - -end - - -""" -Simulation of optimal control given an initial state and an alea scenario - -# Arguments -* `model::SPmodel`: - the DPSPmodel of our problem -* `param::SDPparameters`: - the parameters for the SDP algorithm -* `scenario::Array`: - the scenario of uncertainties realizations we want to simulate on -* `X0::SDPparameters`: - the initial state of the system -* `V::Array`: - the vector representing the value functions as functions of the state - of the system at each time step -* `display::Bool`: - the output display or verbosity parameter - -# Return -* `J::Float`: - the cost of the optimal control over the scenario provided -* `stocks::Array`: - the state of the controlled system at each time step -* `controls::Array`: - the controls applied to the system at each time step -""" -function sdp_forward_single_simulation(model::StochDynProgModel, - param::SDPparameters, - scenario::Array, - X0::Array, - V, - display=true::Bool) - - if VERSION.minor < 5 - scenario = reshape(scenario, size(scenario, 1), size(scenario, 3)) - end - TF = model.stageNumber - law = model.noises - u_bounds = model.ulim - x_bounds = model.xlim - x_steps = param.stateSteps - - #Compute cartesian product spaces - p_states, p_controls = generate_grid(model, param) - product_states = collect(p_states) - product_controls = collect(p_controls) - - product_states = collect(product_states) - product_controls = collect(product_controls) - - controls = Inf*ones(TF-1, 1, model.dimControls) - states = Inf*ones(TF, 1, model.dimStates) - - state_num = 0 - for xj in X0 - state_num += 1 - states[1, 1, state_num] = xj - end - - J = 0 - best_state = X0 - - best_control = tuple() - - if (param.infoStructure == "DH") - #Decision hazard forward simulation - for t = 1:(TF-1) - - x = states[t,1,:] - - best_V = Inf - Vitp = value_function_interpolation(model.dimStates, V, t+1) - - for u in product_controls - - count_admissible_w = 0. - current_V = 0. - - if (param.expectation_computation=="MonteCarlo") - sampling_size = param.monteCarloSize - samples = [sampling(law,t) for i in 1:sampling_size] - probas = (1/sampling_size) - else - sampling_size = law[t].supportSize - samples = law[t].support - probas = law[t].proba - end - - for w = 1:sampling_size - - w_sample = samples[:, w] - proba = probas[w] - - next_state = model.dynamics(t, x, u, w_sample) - - if model.constraints(t, x, u, w_sample)&&SDPutils.is_next_state_feasible(next_state, model.dimStates, model.xlim) - ind_next_state = SDPutils.real_index_from_variable(next_state, x_bounds, x_steps) - next_V = Vitp[ind_next_state...] - current_V += proba *(model.costFunctions(t, x, u, w_sample) + next_V) - count_admissible_w = count_admissible_w + proba - end - end - current_V = current_V/count_admissible_w - if (current_V < best_V)&(count_admissible_w>0) - best_control = u - best_state = model.dynamics(t, x, u, scenario[t,:]) - best_V = current_V - end - end - - index_control = 0 - for uj in best_control - index_control += 1 - controls[t,1,index_control] = uj - end - - index_state = 0 - for xj in best_state - index_state = index_state +1 - states[t+1,1,index_state] = xj - end - J += model.costFunctions(t, x, best_control, scenario[t,:]) - end - - else - #Hazard desision forward simulation - for t = 1:TF-1 - - x = states[t,1,:] - - Vitp = value_function_interpolation(model.dimStates, V, t+1) - - best_V = Inf - - for u = product_controls - - next_state = model.dynamics(t, x, u, scenario[t,:]) - - if model.constraints(t, x, u, scenario[t])&&SDPutils.is_next_state_feasible(next_state, model.dimStates, model.xlim) - ind_next_state = SDPutils.real_index_from_variable(next_state, x_bounds, x_steps) - next_V = Vitp[ind_next_state...] - current_V = model.costFunctions(t, x, u, scenario[t,:]) + next_V - if (current_V < best_V) - best_control = u - best_state = model.dynamics(t, x, u, scenario[t,:]) - best_V = current_V - end - end - - end - index_control = 0 - for uj in best_control - index_control += 1 - controls[t,1,index_control] = uj - end - - index_state = 0 - for xj in best_state - index_state = index_state +1 - states[t+1,1,index_state] = xj - end - J += model.costFunctions(t, x, best_control, scenario[t,:]) - end - end - - x = states[TF,1,:] - J = J + model.finalCostFunction(x) - - return J, states, controls -end - diff --git a/src/SDPutils.jl b/src/SDPutils.jl deleted file mode 100644 index 59cb4ad..0000000 --- a/src/SDPutils.jl +++ /dev/null @@ -1,252 +0,0 @@ -module SDPutils -using Interpolations - -export index_from_variable, real_index_from_variable, - compute_V_given_t_DH, compute_V_given_t_HD - -""" -Convert the state and control float tuples (stored as arrays or tuples) of the -problem into int tuples that can be used as indexes for the discretized -value functions - -# Parameters -* `variable::Array`: - the vector variable we want to convert to an index (integer) -* `bounds::Array`: - the lower bounds for each component of the variable -* `variable_steps::Array`: - discretization step for each component - -# Returns -* `index::Tuple{Integeres}`: - the indexes of the variable - -""" -function index_from_variable(variable, bounds::Array, variable_steps::Array) - return tuple([ 1 + floor(Int64,(1e-10+( variable[i] - bounds[i][1] )/ variable_steps[i] )) for i in 1:length(variable)]...) -end - - -""" -Convert the state and control float tuples (stored as arrays or tuples) of the -problem into float tuples that can be used as indexes for the interpolated -value function - -# Parameters -* `variable::Array`: - the vector variable we want to convert to an index (integer) -* `bounds::Array`: - the lower bounds for each component of the variable -* `variable_steps::Array`: - discretization step for each component - -# Returns -* `index::Tuple{Float64}`: - the indexes of the variable - -""" -function real_index_from_variable(variable, bounds::Array, variable_steps::Array) - return tuple([1 + ( variable[i] - bounds[i][1] )/variable_steps[i] for i in 1:length(variable)]...) -end - -""" -Check if next state x_{t+1} satisfies state bounds constraints - -# Parameters -* `next_stae::Array`: - the state we want to check -* `x_dim::Int`: - the number of state variables -* `x_bounds::Array`: - the state variables bounds - -# Returns -* `index::Tuple{Float64}`: - the indexes of the variable - -""" -function is_next_state_feasible(next_state, x_dim, x_bounds) - - next_state_box_const = true - - for i in 1:x_dim - next_state_box_const = (next_state_box_const&& - (next_state[i]>=x_bounds[i][1])&& - (next_state[i]<=x_bounds[i][2])) - end - - return next_state_box_const -end - -""" -Computes the value function at time t evaluated at state x in a decision -hazard setting - -# Parameters -* `sampling_size::int`: - the size of the uncertainty space -* `samples::Array`: - the uncertainties realizations -* `probas::Array`: - the probabilities of all the uncertainties realizations -* `u_bounds::Array`: - the control variables bounds -* `x_bounds::Array`: - the state variables bounds -* `x_steps::Array`: - the state variables steps -* `x_dim::int`: - the number of state variables -* `product_controls::Array`: - the discretized control space -* `dynamics::Function`: - the dynamic of the problem -* `constraints::Function`: - the constraints of the problem -* `cost::Function`: - the cost function of the problem -* `V::Array`: - the value functions ofthe problem -* `Vitp::Interpolations`: - the interpolated value function at time t+1 -* `t::float`: - the time step at which the value function is computed -* `x::Array or Tuple`: - the state at which the value function needs to be evaluated - -""" -function compute_V_given_x_t_DH(sampling_size, samples, probas, u_bounds, - x_bounds, x_steps, x_dim, product_controls, - dynamics, constraints, cost, V, Vitp, t, x) - expected_V = Inf - optimal_u = tuple() - current_cost = 0 - #Loop over controls - for u = product_controls - - expected_V_u = 0. - count_admissible_w = 0 - - for w = 1:sampling_size - w_sample = samples[:, w] - proba = probas[w] - next_state = dynamics(t, x, u, w_sample) - - if constraints(t, x, u, w_sample)&&is_next_state_feasible(next_state, x_dim, x_bounds) - - count_admissible_w = count_admissible_w + proba - ind_next_state = real_index_from_variable(next_state, x_bounds, - x_steps) - next_V = Vitp[ind_next_state...] - current_cost = cost(t, x, u, w_sample) - expected_V_u += proba*(current_cost + next_V) - - end - end - - if (count_admissible_w>0) - - next_V = next_V / count_admissible_w - - if (expected_V_u < expected_V) - - expected_V = expected_V_u - optimal_u = u - - end - end - end - ind_x = index_from_variable(x, x_bounds, x_steps) - - V[ind_x..., t] = expected_V -end - - -""" -Computes the value function at time t evaluated at state x in a hazard -decision setting - -# Parameters -* `sampling_size::int`: - the size of the uncertainty space -* `samples::Array`: - the uncertainties realizations -* `probas::Array`: - the probabilities of all the uncertainties realizations -* `u_bounds::Array`: - the control variables bounds -* `x_bounds::Array`: - the state variables bounds -* `x_steps::Array`: - the state variables steps -* `x_dim::int`: - the number of state variables -* `product_controls::Array`: - the discretized control space -* `dynamics::Function`: - the dynamic of the problem -* `constraints::Function`: - the constraints of the problem -* `cost::Function`: - the cost function of the problem -* `V::Array`: - the value functions ofthe problem -* `Vitp::Interpolations`: - the interpolated value function at time t+1 -* `t::float`: - the time step at which the value function is computed -* `x::Array or Tuple`: - the state at which the value function needs to be evaluated - -""" -function compute_V_given_x_t_HD(sampling_size, samples, probas, u_bounds, - x_bounds, x_steps, x_dim, product_controls, - dynamics, constraints, cost, V, Vitp, t, x) - - expected_V = 0. - current_cost = 0. - count_admissible_w = 0. - admissible_u_w_count = 0 - best_V_x_w = Inf - next_V_x_w = Inf - - #Compute expectation - for w in 1:sampling_size - admissible_u_w_count = 0 - best_V_x_w = Inf - next_V_x_w = Inf - w_sample = samples[:, w] - proba = probas[w] - - #Loop over controls to find best next value function - for u in product_controls - - next_state = dynamics(t, x, u, w_sample) - - - if constraints(t, x, u, w_sample)&&is_next_state_feasible(next_state, x_dim, x_bounds) - admissible_u_w_count += 1 - current_cost = cost(t, x, u, w_sample) - ind_next_state = real_index_from_variable(next_state, x_bounds, - x_steps) - next_V_x_w_u = Vitp[ind_next_state...] - next_V_x_w = current_cost + next_V_x_w_u - - if (next_V_x_w < best_V_x_w) - best_V_x_w = next_V_x_w - end - - end - end - - expected_V += proba*best_V_x_w - count_admissible_w += (admissible_u_w_count>0)*proba - end - if (count_admissible_w>0.) - expected_V = expected_V / count_admissible_w - end - ind_x = index_from_variable(x, x_bounds, x_steps) - V[ind_x..., t] = expected_V -end - -end diff --git a/src/StochDynamicProgramming.jl b/src/StochDynamicProgramming.jl index da4a16a..e02aace 100644 --- a/src/StochDynamicProgramming.jl +++ b/src/StochDynamicProgramming.jl @@ -1,34 +1,44 @@ -# 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") +include("sdpLoops.jl") using MathProgBase, JuMP, Distributions +using DocStringExtensions +using CutPruners +using Compat export solve_SDDP, NoiseLaw, simulate_scenarios, SDDPparameters, LinearSPModel, set_state_bounds, - extensive_formulation, - PolyhedralFunction, NextStep, forward_simulations, - StochDynProgModel, SDPparameters, solve_DP, - sdp_forward_simulation, sampling, get_control, get_bellman_value, - benchmark_parameters + extensive_formulation, + PolyhedralFunction, forward_simulations, + StochDynProgModel, SDPparameters, solve_dp, + sampling, get_control, get_bellman_value, + benchmark_parameters, SDDPInterface +include("noises.jl") include("objects.jl") +include("stopcrit.jl") +include("params.jl") +include("regularization.jl") +include("interface.jl") include("utils.jl") include("oneStepOneAleaProblem.jl") include("forwardBackwardIterations.jl") include("SDDPoptimize.jl") include("extensiveFormulation.jl") -include("SDPoptimize.jl") +include("sdp.jl") include("compare.jl") include("cutpruning.jl") -include("stoppingtest.jl") +include("simulation.jl") end diff --git a/src/compare.jl b/src/compare.jl index 2db1e85..23a9ac5 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/. @@ -33,7 +33,7 @@ function benchmark_parameters(model, SDDParametersCollection, scenarios::Array{Float64,3}, seeds::Int; - verbose=0) + verbosity=0) #Execute a first time each function to compile them (V, pbs, callsolver), t1, m1 = @timed solve_SDDP(model, SDDParametersCollection[1], 0) @@ -60,19 +60,19 @@ 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) - if verbose > 0 + if verbosity > 0 print("Instance \t") print("Solving time = ",round(solvingtime,4),"\t") print("Solving memory = ", solvingmemory,"\t") 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..930a310 100644 --- a/src/cutpruning.jl +++ b/src/cutpruning.jl @@ -1,41 +1,11 @@ ############################################################################# -# 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 +# Wrapper of cut's pruning algorithm from CutPruners +############################################################################# """ @@ -43,338 +13,48 @@ 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) + 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 - - # 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 +"""Synchronise cuts between `sddp.pruner` and `sddp.bellmanfunctions`.""" +function sync!(sddp::SDDPInterface) + for t in 1:sddp.spmodel.stageNumber-1 + sddp.bellmanfunctions[t] = PolyhedralFunction(sddp.pruner[t].b, sddp.pruner[t].A) 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)) - - -"""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) +"""Prune cuts with exact pruning.""" +function cleancuts!(sddp::SDDPInterface) + ub = [x[2] for x in sddp.spmodel.xlim] + lb = [x[1] for x in sddp.spmodel.xlim] + for t in 1:sddp.spmodel.stageNumber-1 + exactpruning!(sddp.pruner[t], sddp.params.SOLVER, ub=ub, lb=lb) 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 -* `cutscontainer::ActiveCutsContainer`: -* `V::PolyhedralFunction`: - Object storing all cuts -# Return -* `V::PolyhedralFunction` - the new PolyhedralFunction -""" -function level1_cuts_pruning!(model::SPModel, param::SDDPparameters, - V::PolyhedralFunction, cutscontainer::ActiveCutsContainer) - assert(cutscontainer.numCuts == V.numCuts) - - 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)) -end - - -"""Remove useless cuts in PolyhedralFunction - -First test if cuts are active on the visited states, -then test remaining cuts. - -# Arguments -* `model::SPModel`: -* `cutscontainer::ActiveCutsContainer`: -* `V::PolyhedralFunction`: - Object storing all cuts - -# Return -* `V::PolyhedralFunction` - the new PolyhedralFunction -""" -function exact_cuts_pruning_accelerated!(model::SPModel, param::SDDPparameters, - V::PolyhedralFunction, - cutscontainer::ActiveCutsContainer) - - assert(cutscontainer.numCuts == V.numCuts) - solver = param.SOLVER - - nstates = [length(terr) for terr in cutscontainer.territories] - # Set of inactive cuts: - inactive_cuts = nstates .== 0 - # Set of active cuts: - active_cuts = nstates .> 0 - - # get index of inactive cuts: - index = collect(1:cutscontainer.numCuts)[inactive_cuts] - - # Check if inactive cuts are useful or not: - for id in index - status, x = is_cut_relevant(model, id, V, solver) - if status - active_cuts[id] = true - end - end - - # Remove useless cuts: - cutscontainer.territories = cutscontainer.territories[active_cuts] - cutscontainer.numCuts = sum(active_cuts) - return PolyhedralFunction(V.betas[active_cuts], - V.lambdas[active_cuts, :], - sum(active_cuts)) -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]) + +# Update cut pruner +update!{T}(pruner::CutPruners.DeMatosCutPruner, x::Vector{T}, λ::Vector{T})=addposition!(pruner, x) +update!{T}(pruner::CutPruners.AvgCutPruner, x::Vector{T}, λ::Vector{T})=addusage!(pruner, λ) +update!{T}(pruner::CutPruners.DecayCutPruner, x::Vector{T}, λ::Vector{T})=addusage!(pruner, λ) diff --git a/src/extensiveFormulation.jl b/src/extensiveFormulation.jl index abd476a..73c7e9b 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/. @@ -13,14 +13,14 @@ measurability constraints. # Arguments: * `model::SPModel` * `param::SDDPparameters` -* `verbose`::Int` +* `verbosity`::Int` Optionnal, default is 0 # Returns * `objective value` * `first control` * `status of optimization problem` """ -function extensive_formulation(model, param; verbose=0) +function extensive_formulation(model, param; verbosity=0) #Recover all the constant in the model or in param laws = model.noises @@ -103,7 +103,7 @@ function extensive_formulation(model, param; verbose=0) solved = (status == :Optimal) if solved - (verbose > 0) && println("EF value: "*string(getobjectivevalue(mod))) + (verbosity > 0) && println("EF value: "*string(getobjectivevalue(mod))) firstControl = collect(values(getvalue(u)))[1:DIM_CONTROL*laws[1].supportSize] return getobjectivevalue(mod), firstControl, status else diff --git a/src/forwardBackwardIterations.jl b/src/forwardBackwardIterations.jl index 1883628..19a17d9 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,42 @@ 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 + # If regularization 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 + + # run forward pass + costs, stockTrajectories,_,callsolver_forward, tocfw = forward_simulations(model, param, problems_fp, noise_scenarios, - acceleration=param.IS_ACCELERATED) - - model.refTrajectories = stockTrajectories - return costs, stockTrajectories, callsolver_forward + pruner=sddp.pruner, + regularizer=sddp.regularizer, + verbosity = sddp.verbosity) + + # update SDDP's stats + 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 @@ -67,6 +74,12 @@ scenario according to the current value functions. * `xi::Array{float}`: the noise scenarios on which we simulate, each column being one scenario : xi[t,k,:] is the alea at time t of scenario k. +* `pruner::AbstractCutPruner` + Cut pruner +* `regularizer::SDDPRegularization` + SDDP regularization to use in forward pass +* `verbosity::Int` + Log-level # Returns * `costs::Array{float,1}`: @@ -80,15 +93,20 @@ scenario according to the current value functions. scenario k at time t. * `callsolver::Int64`: the number of solver's call' +* `solvertime::Vector{Float64}` + Solver's call execution time """ function forward_simulations(model::SPModel, param::SDDPparameters, solverProblems::Vector{JuMP.Model}, xi::Array{Float64}; - acceleration=false) + pruner=Nullable{AbstractCutPruner}(), + regularizer=Nullable{SDDPRegularization}(), + verbosity::Int64=0) callsolver::Int = 0 + solvertime = Float64[] T = model.stageNumber nb_forward = size(xi)[2] @@ -123,26 +141,40 @@ 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,verbosity = verbosity) else - status, nextstep = solve_one_step_one_alea(model, param, - solverProblems[t], t, state_t, alea_t) + # switch between HD and DH info structure + if model.info == :HD + sol, ts = solve_one_step_one_alea(model, param, + solverProblems[t], t, state_t, alea_t,verbosity=verbosity) + else + sol, ts = solve_dh(model, param, t, state_t, + solverProblems[t],verbosity=verbosity) + end + end + + # update solvertime with ts + push!(solvertime, ts) + + # update cutpruners status with new point + if param.prune && ~isnull(pruner) && t < T-1 + update!(pruner[t+1], sol.xf, sol.πc) end # 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 +185,7 @@ function forward_simulations(model::SPModel, end end end - return costs, stockTrajectories, controls, callsolver + return costs, stockTrajectories, controls, callsolver, solvertime end @@ -161,6 +193,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 +206,17 @@ 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},verbosity=verbosity) + (verbosity > 4) && println("adding cut to polyhedral function at time t=",t) + 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 +isinside(Vt::PolyhedralFunction, lambda::Vector{Float64})=hash(lambda) in Vt.hashcuts """ Add a cut to the JuMP linear problem. @@ -192,21 +230,33 @@ Add a cut to the JuMP linear problem. Time index * `beta::Float`: affine part of the cut to add -* `lambda::Array{float,1}`: +* `lambda::Vector{Float64}`: subgradient of the cut to add """ 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) + t::Int64, beta::Float64, lambda::Vector{Float64},verbosity=verbosity) + (verbosity > 4) && println("adding cut to model at time t=",t) + alpha = problem[:alpha] + xf = problem[:xf] + @constraint(problem, beta + dot(lambda, xf) <= alpha) + problem.ext[:ncuts] += 1 end +function add_cut_dh!(model::SPModel, problem::JuMP.Model, + t::Int64, beta::Float64, lambda::Vector{Float64}, verbosity=verbosity) + (verbosity > 4) && println("adding cut to dh model at time t=",t) + alpha = problem[:alpha] + xf = problem[:xf] + + for j=1:length(model.noises[t].proba) + @constraint(problem, beta + dot(lambda, xf[:, j]) <= alpha[j]) + end +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,92 +264,130 @@ 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}, - stockTrajectories::Array{Float64, 3}, - law) +function backward_pass!(sddp::SDDPInterface, + stockTrajectories::Array{Float64, 3}) - callsolver::Int = 0 + model = sddp.spmodel + law = model.noises + param = sddp.params + solverProblems = sddp.solverinterface + V = sddp.bellmanfunctions + + solvertime = Float64[] T = model.stageNumber nb_forward = size(stockTrajectories)[2] - costs::Vector{Float64} = zeros(1) - state_t = zeros(Float64, model.dimStates) for t = T-1:-1:1 - costs = zeros(Float64, law[t].supportSize) - for k = 1:nb_forward - subgradient_array = zeros(Float64, model.dimStates, law[t].supportSize) # We collect current state: - state_t = collect(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. - proba = zeros(law[t].supportSize) - - # We iterate other the possible realization of noise: - for w in 1:law[t].supportSize - - # We get current noise: - alea_t = collect(law[t].support[:, w]) - - 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) - - if solved - # We catch the subgradient λ: - subgradient_array[:, w] = nextstep.sub_gradient - # and the current cost: - costs[w] = nextstep.cost - # and as problem is solved we store current proba in array: - proba[w] = law[t].proba[w] - end + state_t = stockTrajectories[t, k, :] + if model.info == :HD + compute_cuts_hd!(model, param, V, solverProblems, t, state_t, solvertime,sddp.verbosity) + elseif model.info == :DH + compute_cuts_dh!(model, param, V, solverProblems, t, state_t, solvertime,sddp.verbosity) end - # We add cuts only if one solution was being found: - if sum(proba) > 0 - # Scale probability (useful when some problems where infeasible): - proba /= sum(proba) - - # Compute expectation of subgradient λ: - subgradient = vec(sum(proba' .* subgradient_array, 2)) - # ... expectation of cost: - costs_npass = dot(proba, costs) - # ... and expectation of slope β: - 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) - end + end + end + # update stats + sddp.stats.nsolved += length(solvertime) + sddp.stats.solverexectime_bw = vcat(sddp.stats.solverexectime_bw, solvertime) +end + + +"""Compute cuts in Hazard-Decision (classical SDDP).""" +function compute_cuts_hd!(model::SPModel, param::SDDPparameters, + V::Vector{PolyhedralFunction}, + solverProblems::Vector{JuMP.Model}, t::Int, + state_t::Vector{Float64}, solvertime::Vector{Float64},verbosity::Int64) + law = model.noises + costs = zeros(Float64, model.noises[t].supportSize) + subgradient_array = zeros(Float64, model.dimStates, model.noises[t].supportSize) + + # 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. + proba = zeros(model.noises[t].supportSize) + + # We iterate other the possible realization of noise: + for w in 1:model.noises[t].supportSize + + # We get current noise: + alea_t = collect(model.noises[t].support[:, w]) + # We solve LP problem with current noise and position: + sol, ts = solve_one_step_one_alea(model, param, + solverProblems[t], + t, state_t, alea_t, + relaxation=model.IS_SMIP) + push!(solvertime, ts) + + if sol.status + # We catch the subgradient λ: + subgradient_array[:, w] = sol.ρe + # and the current cost: + costs[w] = sol.objval + # and as problem is solved we store current proba in array: + proba[w] = law[t].proba[w] + end + end + + # We add cuts only if one solution was being found: + if sum(proba) > 0 + # Scale probability (useful when some problems where infeasible): + proba /= sum(proba) + + # Compute expectation of subgradient λ: + subgradient = vec(sum(proba' .* subgradient_array, 2)) + # ... expectation of cost: + costs_npass = dot(proba, costs) + # ... and expectation of slope β: + beta = costs_npass - dot(subgradient, state_t) + + # Add cut to polyhedral function and JuMP model: + if ~isinside(V[t], subgradient) + add_cut!(model, t, V[t], beta, subgradient, verbosity) + if t > 1 + add_cut_to_model!(model, solverProblems[t-1], t, beta, subgradient, verbosity) end + end + end +end + +"""Compute cuts in Decision-Hazard (variant of SDDP).""" +function compute_cuts_dh!(model::SPModel, param::SDDPparameters, + V::Vector{PolyhedralFunction}, + solverProblems::Vector{JuMP.Model}, t::Int, + state_t::Vector{Float64}, solvertime::Vector{Float64},verbosity::Int64) + # We solve LP problem in decision-hazard, considering all possible + # outcomes of randomness: + sol, ts = solve_dh(model, param, t, state_t, solverProblems[t]) + + push!(solvertime, ts) + + # We add cuts only if one solution was being found: + # Scale probability (useful when some problems where infeasible): + if sol.status + # Compute expectation of subgradient λ: + subgradient = sol.ρe + # ... expectation of cost: + costs_npass = sol.objval + # ... and expectation of slope β: + beta = costs_npass - dot(subgradient, state_t) + + # Add cut to polyhedral function and JuMP model: + add_cut!(model, t, V[t], beta, subgradient, verbosity) + if t > 1 + add_cut_dh!(model, solverProblems[t-1], t, beta, subgradient,verbosity) end end - return callsolver end diff --git a/src/interface.jl b/src/interface.jl new file mode 100644 index 0000000..9fe1ea9 --- /dev/null +++ b/src/interface.jl @@ -0,0 +1,79 @@ +# 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} + + verbosity::Int #0: no output, 1: big phases, 2: every verbose_it iterations, 3: inside iterations, 4: detailed inside iterations, 6: showing LP problems + verbose_it::Int + + # Init SDDP interface + function SDDPInterface(model::SPModel, # SP Model + param::SDDPparameters,# parameters + stopcrit::AbstractStoppingCriterion, + prunalgo::AbstractCutPruningAlgo; + regularization=nothing, + verbosity::Int=2, + verbose_it::Int=1) + + check_SDDPparameters(model, param, verbosity) + # initialize value functions: + V, problems = initialize_value_functions(model, param) + (verbosity > 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, verbosity,verbose_it) + end + + function SDDPInterface(model::SPModel, + params::SDDPparameters, + stopcrit::AbstractStoppingCriterion, + prunalgo::AbstractCutPruningAlgo, + V::Vector{PolyhedralFunction}; + regularization=nothing, + verbosity::Int=2, + verbose_it::Int=1) + + check_SDDPparameters(model, params, verbosity) + # 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, verbosity,verbose_it) + end +end + + +function initpruner(algo, n_stages, n_dim) + # Initialize cuts container for cuts pruning: + return [CutPruners.CutPruner{n_dim, Float64}(algo, :Max) for i in 1:n_stages-1] +end + +isregularized(sddp::SDDPInterface) = !isnull(sddp.regularizer) + diff --git a/src/noises.jl b/src/noises.jl index 7a17dbf..7cf33c0 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/. @@ -8,7 +8,6 @@ # - add functions to build scenarios with given probability laws ############################################################################# -using Iterators type NoiseLaw # Number of points in distribution: @@ -87,16 +86,16 @@ function noiselaw_product(law, laws...) nw1 = law.supportSize nw2 = n2.supportSize # and dimensions of aleas: - ndim1 = size(law.support)[1] - ndim2 = size(n2.support)[1] + n_dim1 = size(law.support)[1] + n_dim2 = size(n2.support)[1] # proba and support will defined the output discrete law proba = zeros(nw1*nw2) - support = zeros(ndim1 + ndim2, nw1*nw2) + support = zeros(n_dim1 + n_dim2, nw1*nw2) count = 1 # Use an iterator to find all permutations: - for tup in product(1:nw1, 1:nw2) + for tup in Base.product(1:nw1, 1:nw2) i, j = tup # P(X = (x_i, y_i)) = pi1_i * pi2_i proba[count] = law.proba[i] * n2.proba[j] diff --git a/src/objects.jl b/src/objects.jl index 8f82f9f..6e2d973 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,9 +6,8 @@ # Define all types used in this module. ############################################################################# -include("noises.jl") -abstract SPModel +@compat abstract type SPModel end type PolyhedralFunction @@ -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(n_dim::Int) = PolyhedralFunction(Float64[], Array{Float64}(0, n_dim), 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,15 +51,14 @@ 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} + info::Symbol IS_SMIP::Bool - function LinearSPModel(nstage, # number of stages - ubounds, # bounds of control + function LinearSPModel(n_stage, # number of stages + u_bounds, # bounds of control x0, # initial state cost, # cost function dynamic, # dynamic @@ -59,10 +66,12 @@ type LinearSPModel <: SPModel Vfinal=nothing, # final cost eqconstr=nothing, # equality constraints ineqconstr=nothing, # inequality constraints + info=:HD, # information structure control_cat=nothing) # category of controls + # infer the problem's dimension dimStates = length(x0) - dimControls = length(ubounds) + dimControls = length(u_bounds) dimNoises = length(aleas[1].support[:, 1]) # First step: process terminal costs. @@ -70,26 +79,27 @@ 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] + # control's category + isbu = isa(control_cat, Vector{Symbol})? control_cat: [:Cont for i in 1:dimControls] is_smip = (:Int in isbu)||(:Bin in isbu) - xbounds = [(-Inf, Inf) for i=1:dimStates] + x_bounds = [(-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) + return new(n_stage, dimControls, dimStates, dimNoises, x_bounds, u_bounds, + x0, cost, dynamic, aleas, Vf, isbu, eqconstr, ineqconstr, info, is_smip) end end """Set bounds on state.""" -function set_state_bounds(model::SPModel, xbounds) - if length(xbounds) != model.dimStates +function set_state_bounds(model::SPModel, x_bounds) + if length(x_bounds) != model.dimStates error("Bounds dimension, must be ", model.dimStates) else - model.xlim = xbounds + model.xlim = x_bounds end end @@ -113,6 +123,8 @@ type StochDynProgModel <: SPModel constraints::Function noises::Vector{NoiseLaw} + build_search_space::Nullable{Function} + function StochDynProgModel(model::LinearSPModel, final, cons) if isa(model.costFunctions, Function) cost = model.costFunctions @@ -132,97 +144,14 @@ type StochDynProgModel <: SPModel end function StochDynProgModel(TF, x_bounds, u_bounds, x0, cost_t, - finalCostFunction, dynamic, constraints, aleas) + finalCostFunction, dynamic, constraints, aleas, search_space_builder = Nullable{Function}()) return new(TF, length(u_bounds), length(x_bounds), length(aleas[1].support[:, 1]), x_bounds, u_bounds, x0, cost_t, finalCostFunction, dynamic, - constraints, aleas) + constraints, aleas, search_space_builder) end 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 +191,14 @@ type SDPparameters end -function set_max_iterations(param::SDDPparameters, n_iter::Int) - param.maxItNumber = n_iter -end + +@compat abstract type AbstractSDDPStats 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 +209,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 +245,37 @@ 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 + # cuts' multipliers + πc::Vector{Float64} end +NLDSSolution() = NLDSSolution(false, Inf, Array{Float64, 1}(), Array{Float64, 1}(), Array{Float64, 1}(), Inf, Array{Float64, 1}()) diff --git a/src/oneStepOneAleaProblem.jl b/src/oneStepOneAleaProblem.jl index c2df7b9..15d185c 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/. @@ -38,10 +38,10 @@ problem with respect to the initial state x If specified, approximate future cost as 0 # Returns -* `solved::Bool`: - True if the solution is feasible, false otherwise * `NextStep`: Store solution of the problem +* `ts::Float64`: + Solver's execution time """ function solve_one_step_one_alea(model, param, @@ -50,11 +50,12 @@ function solve_one_step_one_alea(model, xt::Vector{Float64}, xi::Vector{Float64}; relaxation=false::Bool, - init=false::Bool) + init=false::Bool, + verbosity::Int64=0) # Get var defined in JuMP.model: - u = getvariable(m, :u) - w = getvariable(m, :w) - alpha = getvariable(m, :alpha) + u = m[:u] + w = m[:w] + alpha = m[:alpha] # Update value of w: for ii in 1:model.dimNoises @@ -66,61 +67,117 @@ function solve_one_step_one_alea(model, JuMP.setRHS(m.ext[:cons][i], xt[i]) end + if verbosity > 5 + println("One step one alea problem at time t=",t) + println("for x =",xt) + println("and w=",xi) + print(m) + end + if model.IS_SMIP - solved = relaxation ? solve_relaxed!(m, param): solve_mip!(m, param) + solved = relaxation ? solve_relaxed!(m, param,verbosity): solve_mip!(m, param,verbosity) else - status = solve(m, suppress_warnings=true) + status = (verbosity>3) ? solve(m, suppress_warnings=false) : solve(m, suppress_warnings=true) 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), - getvalue(alpha)) + getdual(m.ext[:cons]), + getvalue(alpha), + getcutsmultipliers(m)) + else + # If no solution is found, then return nothing + result = NLDSSolution() + end + + return result, solvetime +end + + +"""Solve model in Decision-Hazard.""" +function solve_dh(model, param, t, xt, m; verbosity::Int64=0) + xf = m[:xf] + u = m[:u] + alpha = m[:alpha] + for i in 1:model.dimStates + JuMP.setRHS(m.ext[:cons][i], xt[i]) + end + + (verbosity>5) && println("Decision Hazard model") + (verbosity>5) && print(m) + + status = solve(m) + solved = status == :Optimal + if ~solved + println(m) + error("Foo") + end + + solvetime = try getsolvetime(m) catch 0 end + + if solved + # Computation of subgradient: + λ = Float64[getdual(m.ext[:cons][i]) for i in 1:model.dimStates] + result = NLDSSolution(solved, + getobjectivevalue(m), + getvalue(xf)[:, 1], + getvalue(u), + λ, + getvalue(alpha)[1], + getcutsmultipliers(m)) else # If no solution is found, then return nothing - result = nothing + result = NLDSSolution() 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},verbosity::Int64=0) # 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)) + xf = m[:xf] + 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) + res = solve_one_step_one_alea(model, param, m, t, xt, xi,verbosity) m.obj = pobj return res end + """Solve relaxed MILP problem.""" -function solve_relaxed!(m, param) +function solve_relaxed!(m, param,verbosity::Int64=0) setsolver(m, param.SOLVER) status = solve(m, relaxation=true) return status == :Optimal end """Solve original MILP problem.""" -function solve_mip!(m, param) - setsolver(m, param.MIPSOLVER) +function solve_mip!(m, param,verbosity::Int64=0) + setsolver(m, get(param.MIPSOLVER)) status = solve(m, relaxation=false) return status == :Optimal end + +getcutsmultipliers(m::JuMP.Model)=_getdual(m)[end-m.ext[:ncuts]+1:end] +function _getdual(m::JuMP.Model) + return MathProgBase.SolverInterface.getconstrduals(m.internalModel) +end diff --git a/src/params.jl b/src/params.jl new file mode 100644 index 0000000..a57580c --- /dev/null +++ b/src/params.jl @@ -0,0 +1,64 @@ +# 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 + # max iterations + max_iterations::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 + # Pruning: + prune::Bool + + 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, prune=false) + + return new(solver, mipsolver, passnumber, max_iterations, confidence, + compute_ub, montecarlo_final, montecarlo_in_iter, reload, prune) + end +end + + +""" +Test compatibility of parameters. + +# Arguments +* `model::SPModel`: + Parametrization of the problem +* `param::SDDPparameters`: + Parameters of SDDP +* `verbosity:Int64`: + +# Return +`Bool` +""" +function check_SDDPparameters(model::SPModel, param::SDDPparameters, verbosity=0::Int64) + if model.IS_SMIP && isnull(param.MIPSOLVER) + error("MIP solver is not defined. Please set `param.MIPSOLVER`") + end + + (verbosity > 0) && (model.IS_SMIP) && println("SMIP SDDP") +end + diff --git a/src/regularization.jl b/src/regularization.jl new file mode 100644 index 0000000..fb0790a --- /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 + +@compat abstract type AbstractRegularization end + + +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/sdp.jl b/src/sdp.jl new file mode 100644 index 0000000..e4e50ff --- /dev/null +++ b/src/sdp.jl @@ -0,0 +1,460 @@ +# 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, Interpolations + + +""" +Compute interpolation of the value function at time t + +# Arguments +* `model::SPmodel`: +* `dim_states::Int`: + the number of state variables +* `v::Array`: + the value function to interpolate +* `time::Int`: + time at which we have to interpolate V + +# Return +* Interpolation + the interpolated value function (working as an array with float indexes) + +""" +function value_function_interpolation( dim_states::Int, V::Union{SharedArray, Array}, time::Int) + return interpolate(V[[Colon() for i in 1:dim_states]...,time], BSpline(Linear()), OnGrid()) +end + + +""" +Compute the cartesian products of discretized state spaces + +# Arguments +* `model::SPmodel`: + the model of the problem +* `param::SDPparameters`: + the parameters of the problem + +# Return +* Iterators: product_states + the cartesian product iterators for states + +""" +function generate_state_grid(model::SPModel, param::SDPparameters, w::Nullable{Array} = Nullable{Array}() ) + product_states = Base.product([model.xlim[i][1]:param.stateSteps[i]:model.xlim[i][2] for i in 1:model.dimStates]...) + + return collect(product_states) +end + +""" +Compute the cartesian products of discretized control spaces or more complex space if provided + +# Arguments +* `model::SPmodel`: + the model of the problem +* `param::SDPparameters`: + the parameters of the problem +* `t::Int`: + time step of the value function computation +* `x::Array{Float64}`: + the current state explored + +# Return +* Iterators: product_states and product_controls + the cartesian product iterators for both states and controls + +""" +function generate_control_grid(model::SPModel, param::SDPparameters, + t::Nullable{Int} = Nullable{Int}(), + x::Nullable{Array} = Nullable{Array}(), + w::Nullable{Array} = Nullable{Array}()) + + if (isnull(model.build_search_space))||(isnull(t))||(isnull(x)) + product_controls = Base.product([model.ulim[i][1]:param.controlSteps[i]:model.ulim[i][2] for i in 1:model.dimControls]...) + else + product_controls = model.build_search_space(t, x, w) + + end + return collect(product_controls) +end + + +""" +Transform a general SPmodel into a StochDynProgModel + +# Arguments +* `model::SPmodel`: + the model of the problem +* `param::SDPparameters`: + the parameters of the problem + +# Return +* `sdpmodel::StochDynProgModel: + the corresponding StochDynProgModel + +""" +function build_sdpmodel_from_spmodel(model::SPModel) + + function zero_fun(x) + return 0 + end + + if isa(model,LinearSPModel) + function cons_fun(t,x,u,w) + return true + end + if in(:finalCostFunction,fieldnames(model)) + SDPmodel = StochDynProgModel(model, model.finalCostFunction, cons_fun) + else + SDPmodel = StochDynProgModel(model, zero_fun, cons_fun) + end + elseif isa(model,StochDynProgModel) + SDPmodel = model + else + error("cannot build StochDynProgModel from current SPmodel. You need to + implement a new StochDynProgModel constructor.") + end + + return SDPmodel +end + + +""" +Dynamic programming algorithm to compute optimal value functions +by backward induction using bellman equation in the finite horizon case. +The information structure can be Decision Hazard (DH) or Hazard Decision (HD) + +# Arguments +* `model::SPmodel`: + the DPSPmodel of our problem +* `param::SDPparameters`: + the parameters for the SDP algorithm +* `display::Int`: + the output display or verbosity parameter + +# Return +* `value_functions::Array`: + the vector representing the value functions as functions of the state + of the system at each time step + +""" +function solve_dp(model::SPModel, param::SDPparameters, display=0::Int64) + + SDPmodel = build_sdpmodel_from_spmodel(model) + + # Start of the algorithm + V = compute_value_functions_grid(SDPmodel, param, display) + return V +end + + +""" +Dynamic Programming algorithm to compute optimal value functions + +# Parameters +* `model::StochDynProgModel`: + the StochDynProgModel of the problem +* `param::SDPparameters`: + the parameters for the algorithm +* `display::Int`: + the output display or verbosity parameter + +# Returns +* `value_functions::Array`: + the vector representing the value functions as functions of the state + of the system at each time step + +""" +function compute_value_functions_grid(model::StochDynProgModel, + param::SDPparameters, + display=0::Int64) + + TF = model.stageNumber + next_state = zeros(Float64, model.dimStates) + + u_bounds = model.ulim + x_bounds = model.xlim + x_steps = param.stateSteps + x_dim = model.dimStates + + dynamics = model.dynamics + constraints = model.constraints + cost = model.costFunctions + + law = model.noises + + build_Ux = model.build_search_space + + #Compute cartesian product spaces + product_states = generate_state_grid(model, param) + + product_controls = generate_control_grid(model, param) + + V = SharedArray{Float64}(zeros(Float64, param.stateVariablesSizes..., TF)) + + #Compute final value functions + for x in product_states + ind_x = SdpLoops.index_from_variable(x, x_bounds, x_steps) + V[ind_x..., TF] = model.finalCostFunction(x) + end + + if param.expectation_computation!="MonteCarlo" && param.expectation_computation!="Exact" + warn("param.expectation_computation should be 'MonteCarlo' or 'Exact'. + Defaulted to 'exact'") + param.expectation_computation="Exact" + end + + if param.infoStructure == "DH" + get_V_t_x = SdpLoops.sdp_u_w_loop + elseif param.infoStructure == "HD" + get_V_t_x = SdpLoops.sdp_w_u_loop + else + warn("Information structure should be DH or HD. Defaulted to DH") + param.infoStructure = "DH" + get_V_t_x = SdpLoops.sdp_u_w_loop + end + + #Construct a progress meter + p = 0 + if display > 0 + p = Progress((TF-1), 1) + println("[SDP] Starting value functions computation:") + end + + # Loop over time: + for t = (TF-1):-1:1 + + if display > 0 + next!(p) + end + + if (param.expectation_computation=="MonteCarlo") + sampling_size = param.monteCarloSize + samples = [sampling(law,t) for i in 1:sampling_size] + probas = (1/sampling_size) + else + sampling_size = law[t].supportSize + samples = law[t].support + probas = law[t].proba + end + + Vitp = value_function_interpolation(x_dim, V, t+1) + + @sync @parallel for indx in 1:length(product_states) + x = product_states[indx] + ind_x = SdpLoops.index_from_variable(x, x_bounds, x_steps) + V[ind_x..., t] = get_V_t_x(sampling_size, samples, probas, + u_bounds, x_bounds, x_steps, x_dim, + product_controls, dynamics, + constraints, cost, Vitp, t, + x, build_Ux)[1] + end + + end + return V +end + +""" +Get the optimal value of the problem from the optimal Bellman Function + +# Arguments +* `model::SPmodel`: + the DPSPmodel of our problem +* `param::SDPparameters`: + the parameters for the SDP algorithm +* `V::Array{Float64}`: + the Bellman Functions + +# Return +* `V_x0::Float64`: + +""" +function get_bellman_value(model::SPModel, param::SDPparameters, + V::Union{SharedArray, Array}) + ind_x0 = SdpLoops.real_index_from_variable(model.initialState, model.xlim, param.stateSteps) + Vi = value_function_interpolation(model.dimStates, V, 1) + println(ind_x0) + return Vi[ind_x0...,1] +end + + +""" +Get the optimal control at time t knowing the state of the system in the decision +hazard case + +# Arguments +* `model::SPmodel`: + the DPSPmodel of our problem +* `param::SDPparameters`: + the parameters for the SDP algorithm +* `V::Array{Float64}`: + the Bellman Functions +* `t::int`: + the time step +* `x::Array`: + the state variable +* `w::Array`: + the alea realization + +# Return +* `V_x0::Float64`: + +""" +function get_control(model::SPModel,param::SDPparameters, + V, t::Int64, x::Array, w::Union{Void, Array} = nothing) + + sdp_model = build_sdpmodel_from_spmodel(model) + + args = [] + optional_args = [] + + if w==nothing + law = sdp_model.noises + get_u = SdpLoops.sdp_dh_get_u + if (param.expectation_computation=="MonteCarlo") + sampling_size = param.monteCarloSize + push!(args,sampling_size, + [sampling(law,t) for i in 1:sampling_size], + (1./sampling_size)*ones(sampling_size)) + else + push!(args,law[t].supportSize, law[t].support, law[t].proba) + end + push!(optional_args, sdp_model.build_search_space) + else + get_u = SdpLoops.sdp_hd_get_u + push!(optional_args, w, sdp_model.build_search_space) + end + + push!(args, sdp_model.ulim, sdp_model.xlim, param.stateSteps, + sdp_model.dimStates, generate_control_grid(sdp_model, param), + sdp_model.dynamics, sdp_model.constraints, sdp_model.costFunctions, + value_function_interpolation(sdp_model.dimStates, V, t+1), t, x) + + return get_u(args..., optional_args...)[1] +end + + +""" +Simulation of optimal control given an initial state and multiple scenarios + +# Arguments +* `model::SPmodel`: + the DPSPmodel of our problem +* `param::SDPparameters`: + the parameters for the SDP algorithm +* `scenarios::Array`: + the scenarios of uncertainties realizations we want to simulate on +* `X0::SDPparameters`: + the initial state of the system +* `V::Array`: + the vector representing the value functions as functions of the state + of the system at each time step +* `display::Bool`: + the output display or verbosity parameter + +# Return +* `costs::Array{Float}`: + the cost of the optimal control over the scenario provided +* `states::Array`: + the state of the controlled system at each time step for each scenario +* `controls::Array`: + the controls applied to the system at each time step for each scenario +""" +function forward_simulations(model::SPModel, + param::SDPparameters, + V::Union{SharedArray, Array}, + scenarios::Array, + display=true::Bool) + + SDPmodel = build_sdpmodel_from_spmodel(model) + + nb_scenarios = size(scenarios,2) + + TF = SDPmodel.stageNumber + law = SDPmodel.noises + x_dim = SDPmodel.dimStates + product_states = generate_state_grid(SDPmodel, param) + costs = SharedArray{Float64}(zeros(nb_scenarios)) + states = SharedArray{Float64}(zeros(TF,nb_scenarios,x_dim)) + controls = SharedArray{Float64}(zeros(TF-1,nb_scenarios,SDPmodel.dimControls)) + + dynamics = SDPmodel.dynamics + cost = SDPmodel.costFunctions + + args = [SDPmodel.ulim, SDPmodel.xlim, param.stateSteps, x_dim, + generate_control_grid(SDPmodel, param), dynamics, SDPmodel.constraints, + cost] + + X0 = SDPmodel.initialState + for s in 1:nb_scenarios + states[1, s, :] = X0 + end + + best_control = tuple() + + info = param.infoStructure + + if info == "DH" + get_u = SdpLoops.sdp_dh_get_u + elseif info == "HD" + get_u = SdpLoops.sdp_hd_get_u + else + warn("Information structure should be DH or HD. Defaulted to DH") + get_u = SdpLoops.sdp_dh_get_u + end + + build_Ux = Nullable{Function}(SDPmodel.build_search_space) + + + @sync @parallel for s in 1:nb_scenarios + + current_scen = scenarios[:,s,:] + + for t = 1:(TF-1) + + args_w = [] + + x = states[t,s,:] + w = current_scen[t,:] + args_t = [value_function_interpolation(x_dim, V, t+1), t, x] + + if info == "DH" + if (param.expectation_computation=="MonteCarlo") + sampling_size = param.monteCarloSize + push!(args_w,sampling_size, + [sampling(law,t) for i in 1:sampling_size], + (1./sampling_size)*ones(sampling_size)) + else + push!(args_w,law[t].supportSize, law[t].support, law[t].proba) + end + else + push!(args_t, w) + end + + best_control = get_u(args_w..., args..., args_t..., build_Ux)[1] + + if best_control == tuple() + error("No u admissible") + else + controls[t,s,:] = [best_control...] + states[t+1,s,:] = dynamics(t, x, best_control, w) + costs[s] = costs[s] + cost(t, x, best_control, w) + end + end + + end + + for s in 1:nb_scenarios + costs[s] = costs[s] + SDPmodel.finalCostFunction(states[TF,s,:]) + end + + return costs, states, controls +end + + diff --git a/src/sdpLoops.jl b/src/sdpLoops.jl new file mode 100644 index 0000000..3659242 --- /dev/null +++ b/src/sdpLoops.jl @@ -0,0 +1,399 @@ +# 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 Bellman equation resolution by +# exhaustive search +# +############################################################################# + +module SdpLoops +using Interpolations + +export index_from_variable, real_index_from_variable + +""" +Convert the state and control float tuples (stored as arrays or tuples) of the +problem into int tuples that can be used as indexes for the discretized +value functions + +# Parameters +* `variable::Array`: + the vector variable we want to convert to an index (integer) +* `bounds::Array`: + the lower bounds for each component of the variable +* `variable_steps::Array`: + discretization step for each component + +# Returns +* `index::Tuple{Integers}`: + the indexes of the variable + +""" +function index_from_variable(variable::Union{Array,Tuple}, bounds::Array, variable_steps::Array) + return tuple([ 1 + floor(Int64,(1e-10+( variable[i] - bounds[i][1] )/ variable_steps[i] )) for i in 1:length(variable)]...) +end + + +""" +Convert the state and control float tuples (stored as arrays or tuples) of the +problem into float tuples that can be used as indexes for the interpolated +value function + +# Parameters +* `variable::Array`: + the vector variable we want to convert to an index (integer) +* `bounds::Array`: + the lower bounds for each component of the variable +* `variable_steps::Array`: + discretization step for each component + +# Returns +* `index::Tuple{Float64}`: + the indexes of the variable + +""" +function real_index_from_variable(variable::Union{Array,Tuple}, bounds::Array, variable_steps::Array) + return tuple([1 + ( variable[i] - bounds[i][1] )/variable_steps[i] for i in 1:length(variable)]...) +end + +""" +Check if next state x_{t+1} satisfies state bounds constraints + +# Parameters +* `next_stae::Array`: + the state we want to check +* `x_dim::Int`: + the number of state variables +* `x_bounds::Array`: + the state variables bounds + +# Returns +* `index::Tuple{Float64}`: + the indexes of the variable + +""" +function is_next_state_feasible(next_state::Union{Array,Tuple}, x_dim::Int, x_bounds::Array) + + next_state_box_const = true + + for i in 1:x_dim + next_state_box_const = (next_state_box_const&& + (next_state[i]>=x_bounds[i][1]-1e-10)&& + (next_state[i]<=x_bounds[i][2]+1e-10)) + end + + return next_state_box_const +end + +""" +Computes the value function at time t evaluated at state x in a decision +hazard setting + +# Parameters +* `sampling_size::int`: + the size of the uncertainty space +* `samples::Array`: + the uncertainties realizations +* `probas::Array`: + the probabilities of all the uncertainties realizations +* `u_bounds::Array`: + the control variables bounds +* `x_bounds::Array`: + the state variables bounds +* `x_steps::Array`: + the state variables steps +* `x_dim::int`: + the number of state variables +* `product_controls::Array`: + the discretized control space +* `dynamics::Function`: + the dynamic of the problem +* `constraints::Function`: + the constraints of the problem +* `cost::Function`: + the cost function of the problem +* `V::Array`: + the value functions ofthe problem +* `Vitp::Interpolations`: + the interpolated value function at time t+1 +* `t::float`: + the time step at which the value function is computed +* `x::Array or Tuple`: + the state at which the value function needs to be evaluated + +# Returns +* `expected_V::Array`: + the value function V(x) +* `optimal_u::Array`: + the optimal control + +""" +function sdp_u_w_loop(sampling_size::Int, samples::Array, + probas::Array, u_bounds::Array, x_bounds::Array, + x_steps::Array, x_dim::Int, product_controls::Array, + dynamics::Function, constraints::Function, cost::Function, + Vitp, t::Int, x::Union{Array,Tuple}, + build_Ux::Nullable{Function} = Nullable{Function}()) + expected_V = Inf + optimal_u = tuple() + #Loop over controls + if isnull(build_Ux) + controls_search_space = product_controls + else + controls_search_space = build_Ux(t,x) + end + + for u in controls_search_space + + expected_V_u = 0. + count_admissible_w = 0 + + for w = 1:sampling_size + w_sample = samples[:, w] + proba = probas[w] + next_state = dynamics(t, x, u, w_sample) + + if constraints(t, x, u, w_sample)&&is_next_state_feasible(next_state, x_dim, x_bounds) + + count_admissible_w = count_admissible_w + proba + + ind_next_state = real_index_from_variable(next_state, x_bounds, + x_steps) + next_V = Vitp[ind_next_state...] + expected_V_u += proba*(cost(t, x, u, w_sample) + next_V) + + end + + end + + if (count_admissible_w>0) + + next_V = next_V / count_admissible_w + + if (expected_V_u < expected_V) + + expected_V = expected_V_u + optimal_u = u + + end + end + end + + return expected_V, optimal_u +end + +""" +Computes the optimal control at time t evaluated +at state x at realization w in a decision hazard setting + +# Parameters +* `sampling_size::int`: + the size of the uncertainty space +* `samples::Array`: + the uncertainties realizations +* `probas::Array`: + the probabilities of all the uncertainties realizations +* `u_bounds::Array`: + the control variables bounds +* `x_bounds::Array`: + the state variables bounds +* `x_steps::Array`: + the state variables steps +* `x_dim::int`: + the number of state variables +* `product_controls::Array`: + the discretized control space +* `dynamics::Function`: + the dynamic of the problem +* `constraints::Function`: + the constraints of the problem +* `cost::Function`: + the cost function of the problem +* `V::Array`: + the value functions ofthe problem +* `Vitp::Interpolations`: + the interpolated value function at time t+1 +* `t::float`: + the time step at which the value function is computed +* `x::Array or Tuple`: + the state at which the value function needs to be evaluated +* `build_Ux::Function or Void`: + an eventual callback to build the controls search at t and x + + +# Returns +* `optimal_u::Array`: + the optimal control + +""" +function sdp_dh_get_u(args...) + + return (sdp_u_w_loop(args...)[2],) + +end + + + +""" +Computes the value function at time t evaluated at state x in a hazard +decision setting + +# Parameters +* `sampling_size::int`: + the size of the uncertainty space +* `samples::Array`: + the uncertainties realizations +* `probas::Array`: + the probabilities of all the uncertainties realizations +* `u_bounds::Array`: + the control variables bounds +* `x_bounds::Array`: + the state variables bounds +* `x_steps::Array`: + the state variables steps +* `x_dim::int`: + the number of state variables +* `product_controls::Array`: + the discretized control space +* `dynamics::Function`: + the dynamic of the problem +* `constraints::Function`: + the constraints of the problem +* `cost::Function`: + the cost function of the problem +* `V::Array`: + the value functions ofthe problem +* `Vitp::Interpolations`: + the interpolated value function at time t+1 +* `t::float`: + the time step at which the value function is computed +* `x::Array or Tuple`: + the state at which the value function needs to be evaluated + + +# Returns +* `expected_V::Array`: + the value function V(x) + +""" +function sdp_w_u_loop(sampling_size::Int, samples::Array, + probas::Array, u_bounds::Array, x_bounds::Array, + x_steps::Array, x_dim::Int, product_controls::Array, + dynamics::Function, constraints::Function, cost::Function, + Vitp, t::Int, x::Union{Array,Tuple}, + build_Ux::Nullable{Function} = Nullable{Function}()) + + expected_V = 0. + count_admissible_w = 0. + + #Compute expectation + for w in 1:sampling_size + admissible_u_w_count = 0 + w_sample = samples[:, w] + proba = probas[w] + + uopt, best_V_x_w, admissible_u_w_count = sdp_hd_get_u(u_bounds, + x_bounds, x_steps, + x_dim, + product_controls, + dynamics, constraints, + cost, Vitp, t, x, w_sample, + build_Ux) + + + expected_V += proba*best_V_x_w + count_admissible_w += (admissible_u_w_count>0)*proba + end + if (count_admissible_w>0.) + expected_V = expected_V / count_admissible_w + end + + return expected_V +end + + +""" +Computes the optimal control and associated cost at time t evaluated +at state x at realization w in a hazard decision setting + +# Parameters +* `u_bounds::Array`: + the control variables bounds +* `x_bounds::Array`: + the state variables bounds +* `x_steps::Array`: + the state variables steps +* `x_dim::int`: + the number of state variables +* `product_controls::Array`: + the discretized control space +* `dynamics::Function`: + the dynamic of the problem +* `constraints::Function`: + the constraints of the problem +* `cost::Function`: + the cost function of the problem +* `V::Array`: + the value functions ofthe problem +* `Vitp::Interpolations`: + the interpolated value function at time t+1 +* `t::float`: + the time step at which the value function is computed +* `x::Array or Tuple`: + the state at which the value function needs to be evaluated +* `w::Array or Tuple`: + the realization at which the value function needs to be evaluated +* `build_Ux::Function or Void`: + an eventual callback to build the controls search at t, x and w + + +# Returns +* `optimal_u::Array`: + the optimal control +* `expected_V::Array`: + the value function V(x) +* `admissible_u_w_count::Int`: + the number of admissible couples (u,w) +""" +function sdp_hd_get_u(u_bounds::Array, x_bounds::Array, + x_steps::Array, x_dim::Int, product_controls::Array, + dynamics::Function, constraints::Function, cost::Function, + Vitp, t::Int, x::Union{Array,Tuple}, w::Union{Array,Tuple}, + build_Ux::Nullable{Function} = Nullable{Function}()) + + if isnull(build_Ux) + controls_search_space = product_controls + else + controls_search_space = get(build_Ux)(t,x,w) + end + + best_V_x_w = Inf + next_V_x_w = Inf + optimal_u = tuple() + admissible_u_w_count = 0 + + for u in product_controls + + next_state = dynamics(t, x, u, w) + + if constraints(t, x, u, w)&&is_next_state_feasible(next_state, x_dim, + x_bounds) + admissible_u_w_count += 1 + ind_next_state = real_index_from_variable(next_state, x_bounds, + x_steps) + next_V_x_w = cost(t, x, u, w) + Vitp[ind_next_state...] + + if (next_V_x_w < best_V_x_w) + best_V_x_w = next_V_x_w + optimal_u = u + end + + end + end + + return optimal_u, best_V_x_w, admissible_u_w_count +end + +end diff --git a/src/stoppingtest.jl b/src/simulation.jl similarity index 72% rename from src/stoppingtest.jl rename to src/simulation.jl index 413da74..cbe1529 100644 --- a/src/stoppingtest.jl +++ b/src/simulation.jl @@ -1,66 +1,42 @@ -# 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 ############################################################################# -""" -Test if the stopping criteria is fulfilled. - -Return true if |upper_bound - lower_bound|/lower_bound < epsilon -or iteration_count > maxItNumber - -# Arguments -*`param::SDDPparameters`: - stopping test type defined in SDDPparameters -* `stats::SDDPStat`: - statistics of the current algorithm - -# Return -`Bool` -""" -function test_stopping_criterion(param::SDDPparameters, stats::SDDPStat) - lb = stats.lower_bounds[end] - ub = stats.upper_bounds[end] + stats.upper_bounds_tol[end] - check_gap = (abs((ub-lb)/lb) < param.gap) - check_iter = stats.niterations >= param.maxItNumber - return check_gap || check_iter -end - - """ Estimate upperbound during SDDP iterations. # Arguments * `model::SPModel`: * `params::SDDPparameters`: -* `Vector{PolyhedralFunction}`: - Polyhedral functions where cuts will be removed * `iteration_count::Int64`: current iteration number * `upperbound_scenarios` -* `verbose::Int64` +* `verbosity::Int64` +* `current_upb::Tuple{Float64}` + Current upper-bound +* `problem::Vector{JuMP.Model}` + Stages' models # Return -* `upb::Float64`: - estimation of upper bound +* `upb, σ, tol`: + estimation of upper bound with confidence level """ function in_iteration_upb_estimation(model::SPModel, param::SDDPparameters, iteration_count::Int64, - verbose::Int64, + verbosity::Int64, upperbound_scenarios, current_upb, problems) upb, σ, tol = current_upb # If specified, compute upper-bound: if (param.compute_ub > 0) && (iteration_count%param.compute_ub==0) - (verbose > 0) && println("Compute upper-bound with ", + (verbosity > 0) && println("Compute upper-bound with ", param.in_iter_mc, " scenarios...") # estimate upper-bound with Monte-Carlo estimation: upb, σ, tol = estimate_upper_bound(model, param, upperbound_scenarios, problems) @@ -85,10 +61,8 @@ Estimate upper bound with Monte Carlo. Number of scenarios to use to compute Monte-Carlo estimation # Return -* `upb::Float64`: +* `upb, σ, tol`: estimation of upper bound -* `costs::Vector{Float64}`: - Costs along different trajectories """ function estimate_upper_bound(model::SPModel, param::SDDPparameters, V::Vector{PolyhedralFunction}, @@ -104,7 +78,7 @@ function estimate_upper_bound(model::SPModel, param::SDDPparameters, problem::Vector{JuMP.Model}) costs = forward_simulations(model, param, problem, aleas)[1] # discard unvalid values: - costs = costs[isfinite(costs)] + costs = costs[isfinite.(costs)] μ = mean(costs) σ = std(costs) tol = upper_bound_confidence(costs, param.confidence_level) @@ -112,6 +86,12 @@ function estimate_upper_bound(model::SPModel, param::SDDPparameters, end +"""Run `nsimu` simulations of SDDP.""" +function simulate(sddp::SDDPInterface, nsimu::Int) + scenarios = simulate_scenarios(sddp.spmodel.noises, nsimu) + forward_simulations(sddp.spmodel, sddp.params, sddp.solverinterface, scenarios)[1:3] +end + """ Estimate the upper bound with a distribution of costs diff --git a/src/stopcrit.jl b/src/stopcrit.jl new file mode 100644 index 0000000..43b3f77 --- /dev/null +++ b/src/stopcrit.jl @@ -0,0 +1,131 @@ +################################################################################ +# SDDP's stopping criterion +# credits to @blegat +# https://github.com/blegat/StochasticDualDynamicProgramming.jl/blob/master/src/stopcrit.jl +################################################################################ +import Base.|, Base.& +export stop, AbstractStoppingCriterion, OrStoppingCriterion, AndStoppingCriterion, IterLimit, Pereira, CutLimit, TimeLimit + +@compat abstract type AbstractStoppingCriterion end + +""" +$(SIGNATURES) + +Returns whether the SDDP algorithm should stop. +If `totalstats.niterations` is 0, no iteration has already been done, otherwise, the `niterations`th iteration has just finished. +This iteration used `stats.npaths` paths and generated `stats.nfcuts` (resp. `stats.nocuts`) new feasibility (resp. optimality) cuts. +The lower bound is now `totalstats.lowerbound` and the upper bound has mean `totalstats.upperbound` and variance `totalstats.σ_UB`. +""" +function stop(s::AbstractStoppingCriterion, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats) + error("`stop' function not defined for $(typeof(s))") +end + +""" +$(TYPEDEF) + +Stops if `lhs` *or* `rhs` want to stop. +""" +type OrStoppingCriterion <: AbstractStoppingCriterion + lhs::AbstractStoppingCriterion + rhs::AbstractStoppingCriterion +end + +function stop(s::OrStoppingCriterion, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats) + stop(s.lhs, stats, totalstats) || stop(s.rhs, stats, totalstats) +end + +function (|)(lhs::AbstractStoppingCriterion, rhs::AbstractStoppingCriterion) + OrStoppingCriterion(lhs, rhs) +end + +""" +$(TYPEDEF) + +Stops if `lhs` *and* `rhs` want to stop. +""" +type AndStoppingCriterion <: AbstractStoppingCriterion + lhs::AbstractStoppingCriterion + rhs::AbstractStoppingCriterion +end + +function stop(s::AndStoppingCriterion, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats) + stop(s.lhs, stats, totalstats) && stop(s.rhs, stats, totalstats) +end + +function (&)(lhs::AbstractStoppingCriterion, rhs::AbstractStoppingCriterion) + AndStoppingCriterion(lhs, rhs) +end + +""" +$(TYPEDEF) + +Stops if `iter` ≧ `limit`. +""" +type IterLimit <: AbstractStoppingCriterion + limit::Int +end + +function stop(s::IterLimit, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats) + totalstats.niterations >= s.limit +end + +""" +$(TYPEDEF) + +Stops if there was less than or equal to `limit` cuts added in the iteration. +For instance, `CutLimit(0)` stops when there are no cuts added. +""" +type CutLimit <: AbstractStoppingCriterion + limit::Int +end + +function stop(s::CutLimit, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats) + totalstats.niterations > 0 && stats.nfcuts + stats.nocuts <= s.limit +end + + +""" +$(TYPEDEF) + +Stops if total time of execution is greater than the time limit specified. +For instance, `TimeLimit(100)` stops after 100s. +""" +type TimeLimit <: AbstractStoppingCriterion + timelimit::Float64 +end + +function stop(s::TimeLimit, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats) + totalstats.niterations > 0 && totalstats.time > s.timelimit +end + + +""" +$(TYPEDEF) + +Stops if `z_UB - α * σ/√K - tol < z_LB < z_UB + α * σ/√K + tol` and `σ / √K > β * max(1, |z_LB|))` +""" +type Pereira <: AbstractStoppingCriterion + α::Float64 + β::Float64 + tol::Float64 + + Pereira(α=2.0, β=0.05, tol=1e-6) = new(α, β, tol) +end + +function stop(s::Pereira, stats::AbstractSDDPStats, totalstats::AbstractSDDPStats) + z_UB = stats.upperbound + z_LB = stats.lowerbound + K = stats.npaths + σ = stats.σ_UB + + if totalstats.niterations > 0 + @assert K >= 0 + σ1 = σ / √K + # On the test optimize_stock with Clp, z_LB = -2, z_UB = -1.999999999999 and σ1 = 0 + # this shows the necessicity for a tolerance + σ2 = s.α * σ1 + s.tol + z_UB - σ2 <= z_LB <= z_UB + σ2 && σ1 < s.β * max(1, abs(z_LB)) + else + false + end +end diff --git a/src/utils.jl b/src/utils.jl index 3fffd2b..d171d96 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,55 +67,13 @@ 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 end -"""Concatenate collection of arrays of PolyhedralFunction.""" -function catcutsarray(polyfunarray::Vector{StochDynamicProgramming.PolyhedralFunction}...) - assert(length(polyfunarray) > 0) - ntimes = length(polyfunarray[1]) - # Concatenate cuts in polyfunarray, and discard final time as we do not add cuts at final time: - concatcuts = StochDynamicProgramming.PolyhedralFunction[catcuts([V[t] for V in polyfunarray]...) for t in 1:ntimes-1] - return vcat(concatcuts, polyfunarray[1][end]) -end - - -"""Concatenate collection of PolyhedralFunction.""" -function catcuts(Vts::StochDynamicProgramming.PolyhedralFunction...) - betas = vcat([V.betas for V in Vts]...) - lambdas = vcat([V.lambdas for V in Vts]...) - numcuts = sum([V.numCuts for V in Vts]) - 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 +90,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..86e7bda 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, CutPruners +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,145 @@ 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] - 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 + @testset "Decision-Hazard" begin + model_dh = StochDynamicProgramming.LinearSPModel(n_stages, u_bounds, + x0, cost, dynamic, laws, + info=:DH) + set_state_bounds(model_dh, x_bounds) + param_dh = StochDynamicProgramming.SDDPparameters(solver, + passnumber=1, + gap=0.001, + max_iterations=10) + sddp_dh = solve_SDDP(model_dh, param_dh, 1) + end + + @testset "Cut-pruning" begin + param_pr = StochDynamicProgramming.SDDPparameters(solver, + passnumber=1, + gap=0.001, + reload=2, prune=true) + + sddppr = SDDPInterface(model, param_pr, + StochDynamicProgramming.IterLimit(10), + CutPruners.DeMatosPruningAlgo(-1), + verbosity=0) + # solve SDDP + solve!(sddppr) + + # test exact cuts pruning + ncutini = StochDynamicProgramming.ncuts(sddppr.bellmanfunctions) + StochDynamicProgramming.cleancuts!(sddppr) + @test StochDynamicProgramming.ncuts(sddppr.bellmanfunctions) <= ncutini end - context("Quadratic regularization") do + @testset "Quadratic regularization" begin 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) + alpha = 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 +191,46 @@ 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 +273,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 +309,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..5ebe2dd 100644 --- a/test/sdp.jl +++ b/test/sdp.jl @@ -1,34 +1,34 @@ ################################################################################ # Test SDDP functions ################################################################################ -using FactCheck, StochDynamicProgramming -using StochDynamicProgramming.SDPutils +using Base.Test, StochDynamicProgramming +using StochDynamicProgramming.SdpLoops -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] var = [0.4, 3.7, 1.9] vart = [0.42, 3.78, 1.932] - ind = SDPutils.index_from_variable(var, bounds, steps) - ind2 = SDPutils.real_index_from_variable(vart, bounds, steps) + ind = SdpLoops.index_from_variable(var, bounds, steps) + ind2 = SdpLoops.real_index_from_variable(vart, bounds, steps) - checkFalse = SDPutils.is_next_state_feasible([0,1,2],3,bounds) - checkTrue = SDPutils.is_next_state_feasible([0.12,1.3,1.3],3,bounds) + checkFalse = SdpLoops.is_next_state_feasible([0,1,2],3,bounds) + checkTrue = SdpLoops.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 @@ -112,8 +112,8 @@ facts("SDP algorithm") do aleas_scen = zeros(2, 1, 1) aleas_scen[:, 1, 1] = alea_year; - stateSteps = [1,1]; - controlSteps = [1,1]; + stateSteps = [2,2]; + controlSteps = [2,2]; monteCarloSize = 2; modelSDP = StochDynProgModel(TF, x_bounds, u_bounds, @@ -123,11 +123,11 @@ facts("SDP algorithm") do paramsSDP = StochDynamicProgramming.SDPparameters(modelSDP, stateSteps, controlSteps, - infoStruct, + "HD", "Exact"); - context("Compare StochDynProgModel constructors") do + @testset "Compare StochDynProgModel constructors" begin modelSDPPiecewise = StochDynamicProgramming.LinearSPModel(TF, @@ -160,38 +160,32 @@ 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); + solve_dp(modelSDP, paramsSDP, false); + @test paramsSDP.infoStructure == "DH" paramsSDP.infoStructure = infoStruct - V_sdp = solve_DP(modelSDP, paramsSDP, false); + 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, - aleas_scen, x0, - V_sdp, true ) - - costs_sdp2, stocks_sdp2, controls_sdp2 = StochDynamicProgramming.sdp_forward_simulation(modelSDP, + costs_sdp2, stocks_sdp2, controls_sdp2 = StochDynamicProgramming.forward_simulations(modelSDP, paramsSDP, - aleas_scen, - V_sdp, true ) - - @fact costs_sdp2[1] --> costs_sdp + V_sdp, + aleas_scen) x = x0 - V_sdp2 = StochDynamicProgramming.sdp_compute_value_functions(modelSDP, paramsSDP, false); + V_sdp2 = StochDynamicProgramming.compute_value_functions_grid(modelSDP, paramsSDP, false); paramsSDP.infoStructure = "DH" - V_sdp3 = StochDynamicProgramming.sdp_compute_value_functions(modelSDP, paramsSDP, false); + V_sdp3 = StochDynamicProgramming.compute_value_functions_grid(modelSDP, paramsSDP, false); paramsSDP.infoStructure = "HD" Vitp = StochDynamicProgramming.value_function_interpolation( modelSDP.dimStates, V_sdp, 1) @@ -202,19 +196,20 @@ 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, + costs_sdp3, stocks_sdp3, controls_sdp3 = StochDynamicProgramming.forward_simulations(modelSDP, paramsSDP, - aleas_scen, - V_sdp3, true ) + V_sdp3, + aleas_scen) 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) + a = StochDynamicProgramming.generate_state_grid(modelSDP, paramsSDP) + b = StochDynamicProgramming.generate_control_grid(modelSDP, paramsSDP) x_bounds = modelSDP.xlim x_steps = paramsSDP.stateSteps @@ -222,33 +217,31 @@ 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]) - - ind = SDPutils.index_from_variable(x, x_bounds, x_steps) - @fact get_bellman_value(modelSDP, paramsSDP, V_sdp2) --> V_sdp2[ind...,1] + @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]) - @fact size(V_sdp) --> (paramsSDP.stateVariablesSizes..., TF) - @fact V_sdp2[1,1,1] <= V_sdp3[1,1,1] --> true + modelSDP.initialState = [xi[1] for xi in x_bounds] + ind = SdpLoops.index_from_variable(modelSDP.initialState, x_bounds, x_steps) + @test get_bellman_value(modelSDP, paramsSDP, V_sdp2) == V_sdp2[ind...,1] + modelSDP.initialState = x0 - @fact size(stocks_sdp) --> (3,1,2) - @fact size(controls_sdp) --> (2,1,2) + @test size(V_sdp) == (paramsSDP.stateVariablesSizes..., TF) + @test V_sdp2[1,1,1] <= V_sdp3[1,1,1] state_ref = zeros(2) - state_ref[1] = stocks_sdp[2,1,1] - state_ref[2] = stocks_sdp[2,1,2] + state_ref[1] = stocks_sdp2[2,1,1] + state_ref[2] = stocks_sdp2[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 (get_control(modelSDP,paramsSDP,V_sdp3, 1, x, w)[1] >= CONTROL_MIN) == true + @test (get_control(modelSDP,paramsSDP,V_sdp3, 1, x, w)[1] <= CONTROL_MAX) == true 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_sdp2) == (3,1,2) + @test size(controls_sdp2) == (2,1,2) end