|
| 1 | +# Copyright 2015, Vincent Leclere, Francois Pacaud and Henri Gerard |
| 2 | +# This Source Code Form is subject to the terms of the Mozilla Public |
| 3 | +# License, v. 2.0. If a copy of the MPL was not distributed with this |
| 4 | +# file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| 5 | +############################################################################# |
| 6 | +# Test SDDP with dam example |
| 7 | +# Source: Adrien Cassegrain |
| 8 | +############################################################################# |
| 9 | + |
| 10 | +#srand(2713) |
| 11 | +push!(LOAD_PATH, "../src") |
| 12 | + |
| 13 | +using StochDynamicProgramming, JuMP, Clp |
| 14 | + |
| 15 | +#Constant that the user have to define himself |
| 16 | +const SOLVER = ClpSolver() |
| 17 | +# const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0) |
| 18 | + |
| 19 | +const N_STAGES = 5 |
| 20 | +const N_SCENARIOS = 3 |
| 21 | + |
| 22 | +const DIM_STATES = 2 |
| 23 | +const DIM_CONTROLS = 2 |
| 24 | +const DIM_ALEAS = 1 |
| 25 | + |
| 26 | + |
| 27 | + |
| 28 | +#Constants that the user does not have to define |
| 29 | +const T0 = 1 |
| 30 | + |
| 31 | +# Constants: |
| 32 | +const VOLUME_MAX = 100 |
| 33 | +const VOLUME_MIN = 0 |
| 34 | + |
| 35 | +const CONTROL_MAX = round(Int, .4/7. * VOLUME_MAX) + 1 |
| 36 | +const CONTROL_MIN = 0 |
| 37 | + |
| 38 | +const W_MAX = round(Int, .5/7. * VOLUME_MAX) |
| 39 | +const W_MIN = 0 |
| 40 | +const DW = 1 |
| 41 | + |
| 42 | +# Define aleas' space: |
| 43 | +const N_ALEAS = Int(round(Int, (W_MAX - W_MIN) / DW + 1)) |
| 44 | +const ALEAS = linspace(W_MIN, W_MAX, N_ALEAS) |
| 45 | + |
| 46 | +const EPSILON = .05 |
| 47 | +const MAX_ITER = 20 |
| 48 | + |
| 49 | +alea_year = Array([7.0 7.0 8.0 3.0 1.0 1.0 3.0 4.0 3.0 2.0 6.0 5.0 2.0 6.0 4.0 7.0 3.0 4.0 1.0 1.0 6.0 2.0 2.0 8.0 3.0 7.0 3.0 1.0 4.0 2.0 4.0 1.0 3.0 2.0 8.0 1.0 5.0 5.0 2.0 1.0 6.0 7.0 5.0 1.0 7.0 7.0 7.0 4.0 3.0 2.0 8.0 7.0]) |
| 50 | + |
| 51 | +const X0 = [50, 50] |
| 52 | + |
| 53 | +Ax=[] |
| 54 | +Au=[] |
| 55 | +Aw=[] |
| 56 | + |
| 57 | +Cx=[] |
| 58 | +Cu=[] |
| 59 | +Cw=[] |
| 60 | + |
| 61 | +function generate_random_dynamic() |
| 62 | + for i=1:N_STAGES |
| 63 | + push!(Ax, rand(DIM_STATES,DIM_STATES)) |
| 64 | + push!(Au, rand(DIM_STATES,DIM_CONTROLS)) |
| 65 | + push!(Aw, rand(DIM_STATES,DIM_ALEAS)) |
| 66 | + end |
| 67 | +end |
| 68 | + |
| 69 | +function generate_random_costs() |
| 70 | + for i=1:N_STAGES |
| 71 | + push!(Cx, rand(1,DIM_STATES)) |
| 72 | + push!(Cu, -1*rand(1,DIM_CONTROLS)) |
| 73 | + push!(Cw, rand(1,DIM_ALEAS)) |
| 74 | + end |
| 75 | +end |
| 76 | + |
| 77 | +generate_random_dynamic() |
| 78 | +generate_random_costs() |
| 79 | + |
| 80 | +# Define dynamic of the dam: |
| 81 | +function dynamic(t, x, u, w) |
| 82 | + return Ax[t]*x+Au[t]*u+Aw[t]*w |
| 83 | +end |
| 84 | + |
| 85 | +# Define cost corresponding to each timestep: |
| 86 | +function cost_t(t, x, u, w) |
| 87 | + return (Cx[t]*x)[1,1]+(Cu[t]*u)[1,1]+(Cw[t]*w)[1,1] |
| 88 | +end |
| 89 | + |
| 90 | + |
| 91 | +"""Build aleas probabilities for each month.""" |
| 92 | +function build_aleas() |
| 93 | + aleas = zeros(N_ALEAS, N_STAGES) |
| 94 | + |
| 95 | + # take into account seasonality effects: |
| 96 | + unorm_prob = linspace(1, N_ALEAS, N_ALEAS) |
| 97 | + proba1 = unorm_prob / sum(unorm_prob) |
| 98 | + proba2 = proba1[N_ALEAS:-1:1] |
| 99 | + |
| 100 | + for t in 1:N_STAGES |
| 101 | + aleas[:, t] = (1 - sin(pi*t/N_STAGES)) * proba1 + sin(pi*t/N_STAGES) * proba2 |
| 102 | + end |
| 103 | + return aleas |
| 104 | +end |
| 105 | + |
| 106 | + |
| 107 | +"""Build an admissible scenario for water inflow.""" |
| 108 | +function build_scenarios(n_scenarios::Int64, probabilities) |
| 109 | + scenarios = zeros(n_scenarios, N_STAGES) |
| 110 | + |
| 111 | + for scen in 1:n_scenarios |
| 112 | + for t in 1:N_STAGES |
| 113 | + Pcum = cumsum(probabilities[:, t]) |
| 114 | + |
| 115 | + n_random = rand() |
| 116 | + prob = findfirst(x -> x > n_random, Pcum) |
| 117 | + scenarios[scen, t] = prob |
| 118 | + end |
| 119 | + end |
| 120 | + return scenarios |
| 121 | +end |
| 122 | + |
| 123 | + |
| 124 | +"""Build probability distribution at each timestep. |
| 125 | +
|
| 126 | +Return a Vector{NoiseLaw}""" |
| 127 | +function generate_probability_laws() |
| 128 | + aleas = build_scenarios(N_SCENARIOS, build_aleas()) |
| 129 | + |
| 130 | + laws = Vector{NoiseLaw}(N_STAGES) |
| 131 | + |
| 132 | + # uniform probabilities: |
| 133 | + proba = 1/N_SCENARIOS*ones(N_SCENARIOS) |
| 134 | + |
| 135 | + for t=1:N_STAGES |
| 136 | + laws[t] = NoiseLaw(aleas[:, t], proba) |
| 137 | + end |
| 138 | + |
| 139 | + return laws |
| 140 | +end |
| 141 | + |
| 142 | + |
| 143 | +"""Instantiate the problem.""" |
| 144 | +function init_problem() |
| 145 | + |
| 146 | + x0 = X0 |
| 147 | + aleas = generate_probability_laws() |
| 148 | + |
| 149 | + x_bounds = [(VOLUME_MIN, VOLUME_MAX), (VOLUME_MIN, VOLUME_MAX)] |
| 150 | + u_bounds = [] |
| 151 | + for u = 1:DIM_CONTROLS |
| 152 | + u_bounds = push!(u_bounds, (CONTROL_MIN, CONTROL_MAX)) |
| 153 | + end |
| 154 | + |
| 155 | + model = LinearDynamicLinearCostSPmodel(N_STAGES, |
| 156 | + u_bounds, |
| 157 | + x0, |
| 158 | + cost_t, |
| 159 | + dynamic, |
| 160 | + aleas) |
| 161 | + |
| 162 | + #set_state_bounds(model, x_bounds) |
| 163 | + |
| 164 | + solver = SOLVER |
| 165 | + params = SDDPparameters(solver, N_SCENARIOS, EPSILON, MAX_ITER) |
| 166 | + |
| 167 | + return model, params |
| 168 | +end |
| 169 | + |
| 170 | +model, params = init_problem() |
| 171 | +modelbis = deepcopy(model) |
| 172 | +paramsbis = deepcopy(params) |
| 173 | + |
| 174 | + |
| 175 | +"""Solve the problem.""" |
| 176 | +function solve_dams(model,params,display=false) |
| 177 | + |
| 178 | + V, pbs = solve_SDDP(model, params, display) |
| 179 | + |
| 180 | + aleas = simulate_scenarios(model.noises, |
| 181 | + (model.stageNumber, |
| 182 | + params.forwardPassNumber, |
| 183 | + model.dimNoises)) |
| 184 | + |
| 185 | + params.forwardPassNumber = 1 |
| 186 | + |
| 187 | + costs, stocks = forward_simulations(model, params, V, pbs, aleas) |
| 188 | + println("SDDP cost: ", costs) |
| 189 | + return stocks, V |
| 190 | +end |
| 191 | + |
| 192 | +#Solve the problem and try nb_iter times to generate radom data in case of infeasibility |
| 193 | +unsolve = true |
| 194 | +sol = 0 |
| 195 | +i = 0 |
| 196 | +nb_iter = 10 |
| 197 | + |
| 198 | +while i<nb_iter |
| 199 | + sol, status = extensive_formulation(model,params) |
| 200 | + if (status == :Optimal) |
| 201 | + unsolve = false |
| 202 | + break |
| 203 | + else |
| 204 | + println("\nGenerate new dynamic to reach feasability\n") |
| 205 | + Ax=[] |
| 206 | + Au=[] |
| 207 | + Aw=[] |
| 208 | + generate_random_dynamic() |
| 209 | + model, params = init_problem() |
| 210 | + modelbis = model |
| 211 | + paramsbis = params |
| 212 | + i = i+1 |
| 213 | + end |
| 214 | +end |
| 215 | + |
| 216 | + |
| 217 | +if (unsolve) |
| 218 | + println("Change your parameters") |
| 219 | +else |
| 220 | + a,b = solve_dams(modelbis,paramsbis) |
| 221 | + println("solution =",sol) |
| 222 | +end |
| 223 | + |
0 commit comments