Skip to content

Commit 134acb9

Browse files
frapacfrapac
authored
Merge pull request #175 from JuliaStochOpt/fp/fix_examples
Fix a bug in SDDP after update to JuMP 0.19
2 parents f667f51 + 5edc7b1 commit 134acb9

17 files changed

+160
-224
lines changed

‎Project.toml‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -4,23 +4,23 @@ authors = ["mrchaos <[email protected]>"]
44
version = "0.5.0"
55

66
[deps]
7+
Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
8+
Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
79
CutPruners = "65d46eb8-70e9-5a30-bf48-2afa3a021b8f"
10+
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
11+
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
812
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
913
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
1014
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
11-
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
1215
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
13-
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
14-
Nullables = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
15-
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
1616
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
17-
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
18-
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
17+
Nullables = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
1918
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
19+
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
20+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
2021
SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
21-
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
22-
Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
23-
Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
22+
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
23+
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
2424

2525
[compat]
2626
julia = "^1.3.0"

‎examples/battery_storage_parallel.jl‎

Lines changed: 1 addition & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -29,14 +29,10 @@
2929
# (d) : we don't sell electricity
3030

3131

32-
33-
3432
using StochDynamicProgramming, Distributions
3533
using Statistics
3634
using Distributed
3735

38-
println("library loaded")
39-
4036
# We have to define the instance on all the workers (processes)
4137
@everywhere begin
4238

@@ -104,14 +100,13 @@ println("library loaded")
104100
stateSteps = [0.01]
105101
controlSteps = [0.001]
106102
infoStruct = "HD" # noise at time t is not known before taking the decision at time t
107-
108103
paramSDP = StochDynamicProgramming.SDPparameters(spmodel, stateSteps,
109104
controlSteps, infoStruct)
110105
end
111106

112107
Vs = StochDynamicProgramming.solve_dp(spmodel,paramSDP, 1)
113108

114109
lb_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs)
115-
println("Value obtained by SDP: "*string(lb_sdp))
110+
println("Value obtained by SDP: ", lb_sdp)
116111
costsdp, states, stocks = StochDynamicProgramming.forward_simulations(spmodel,paramSDP,Vs,scenarios)
117112
println(mean(costsdp))

‎examples/dam.jl‎

Lines changed: 10 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -7,10 +7,10 @@
77
# Source: Adrien Cassegrain
88
#############################################################################
99

10-
1110
using StochDynamicProgramming, JuMP
1211

1312
include("solver.jl")
13+
1414
const EPSILON = .05
1515
const MAX_ITER = 20
1616

@@ -20,35 +20,26 @@ const N_SCENARIOS = 10
2020

2121
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])
2222

23-
2423
# COST:
2524
const COST = -66*2.7*(1 .+ .5*(rand(N_STAGES) .- .5))
26-
2725
# Constants:
2826
const VOLUME_MAX = 100
2927
const VOLUME_MIN = 0
30-
3128
const CONTROL_MAX = round(Int, .4/7. * VOLUME_MAX) + 1
3229
const CONTROL_MIN = 0
33-
3430
const W_MAX = round(Int, .5/7. * VOLUME_MAX)
3531
const W_MIN = 0
3632
const DW = 1
37-
3833
const T0 = 1
3934
const HORIZON = 52
40-
4135
const X0 = [90]
4236

4337
# Define aleas' space:
4438
const N_ALEAS = Int(round(Int, (W_MAX - W_MIN) / DW + 1))
4539
const ALEAS = range(W_MIN,stop=W_MAX,length=N_ALEAS)
4640

47-
4841
# Define dynamic of the dam:
49-
function dynamic(t, x, u, w)
50-
return [x[1] - u[1] - u[2] + w[1]]
51-
end
42+
dynamic(t, x, u, w) = [x[1] - u[1] - u[2] + w[1]]
5243

5344
# Define cost corresponding to each timestep:
5445
function cost_t(t, x, u, w)
@@ -59,23 +50,16 @@ end
5950
"""Solve the problem assuming the aleas are known
6051
in advance."""
6152
function solve_determinist_problem()
62-
println(alea_year)
6353
m = JuMP.Model(OPTIMIZER)
64-
6554
@variable(m, 0 <= x[1:N_STAGES] <= 100)
6655
@variable(m, 0. <= u[1:N_STAGES-1] <= 7)
6756
@variable(m, 0. <= s[1:N_STAGES-1] <= 7)
68-
6957
@objective(m, Min, sum(COST[i]*u[i] for i = 1:N_STAGES-1))
70-
7158
for i in 1:(N_STAGES-1)
7259
@constraint(m, x[i+1] - x[i] + u[i] + s[i] - alea_year[i] == 0)
7360
end
74-
7561
@constraint(m, x[1] .==X0)
76-
7762
status = JuMP.optimize!(m)
78-
println(status)
7963
println(JuMP.getobjectivevalue(m))
8064
return JuMP.value.(u), JuMP.value.(x)
8165
end
@@ -84,12 +68,10 @@ end
8468
"""Build aleas probabilities for each month."""
8569
function build_aleas()
8670
aleas = zeros(N_ALEAS, N_STAGES)
87-
8871
# take into account seasonality effects:
8972
unorm_prob = range(1,stop=N_ALEAS,length=N_ALEAS)
9073
proba1 = unorm_prob / sum(unorm_prob)
9174
proba2 = proba1[N_ALEAS:-1:1]
92-
9375
for t in 1:(N_STAGES-1)
9476
aleas[:, t] = (1 - sin(pi*t/N_STAGES)) * proba1 + sin(pi*t/N_STAGES) * proba2
9577
end
@@ -104,7 +86,6 @@ function build_scenarios(n_scenarios::Int64, probabilities)
10486
for scen in 1:n_scenarios
10587
for t in 1:(N_STAGES-1)
10688
Pcum = cumsum(probabilities[:, t])
107-
10889
n_random = rand()
10990
prob = findfirst(x -> x > n_random, Pcum)
11091
scenarios[scen, t] = prob
@@ -119,16 +100,12 @@ end
119100
Return a Vector{NoiseLaw}"""
120101
function generate_probability_laws()
121102
aleas = build_scenarios(N_SCENARIOS, build_aleas())
122-
123-
laws = Vector{NoiseLaw}(N_STAGES-1)
124-
103+
laws = NoiseLaw[]
125104
# uniform probabilities:
126105
proba = 1/N_SCENARIOS*ones(N_SCENARIOS)
127-
128106
for t=1:(N_STAGES-1)
129-
laws[t] = NoiseLaw(aleas[:, t], proba)
107+
push!(laws, NoiseLaw(aleas[:, t], proba))
130108
end
131-
132109
return laws
133110
end
134111

@@ -138,35 +115,24 @@ function init_problem()
138115
# Instantiate model:
139116
x0 = X0
140117
aleas = generate_probability_laws()
141-
142118
x_bounds = [(0, 100)]
143119
u_bounds = [(0, 7), (0, 7)]
144120
model = StochDynamicProgramming.LinearSPModel(N_STAGES,
145121
u_bounds,
146122
x0,
147123
cost_t,
148124
dynamic, aleas)
149-
150125
set_state_bounds(model, x_bounds)
151-
152126
optimizer = OPTIMIZER
153127
params = StochDynamicProgramming.SDDPparameters(optimizer,
154-
passnumber=N_SCENARIOS,
128+
passnumber=1,
155129
gap=EPSILON,
156130
max_iterations=MAX_ITER)
157-
158131
return model, params
159132
end
160133

161-
162-
"""Solve the problem."""
163-
function solve_dams(display=0)
164-
model, params = init_problem()
165-
166-
sddp = solve_SDDP(model, params, display)
167-
aleas = simulate_scenarios(model.noises, params.forwardPassNumber)
168-
169-
costs, stocks = forward_simulations(model, params, sddp.solverinterface, aleas)
170-
println("SDDP cost: ", costs)
171-
return stocks
172-
end
134+
model, params = init_problem()
135+
sddp = solve_SDDP(model, params, 1)
136+
aleas = simulate_scenarios(model.noises, params.forwardPassNumber)
137+
costs, stocks = forward_simulations(model, params, sddp.solverinterface, aleas)
138+
println("SDDP cost: ", costs)

‎examples/damsvalley.jl‎

Lines changed: 21 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -9,28 +9,21 @@ using StochDynamicProgramming, JuMP, LinearAlgebra
99
Random.seed!(2713)
1010

1111
include("solver.jl")
12-
##################################################
13-
1412

1513
##################################################
1614
# PROBLEM DEFINITION
1715
##################################################
1816
# We consider here a valley with 5 dams:
1917
const N_DAMS = 5
20-
2118
const N_STAGES = 12
2219
const N_ALEAS = 10
23-
2420
# Cost are negative as we sell the electricity produced by
2521
# dams (and we want to minimize our problem)
2622
const COST = -66*2.7*(1 .+ .5*(rand(N_STAGES) .- .5))
27-
28-
# Constants:
29-
const VOLUME_MAX = 80
30-
const VOLUME_MIN = 0
31-
32-
const CONTROL_MAX = 40
33-
const CONTROL_MIN = 0
23+
const VOLUME_MAX = 80.0
24+
const VOLUME_MIN = 0.0
25+
const CONTROL_MAX = 40.0
26+
const CONTROL_MIN = 0.0
3427

3528
# Define initial status of stocks:
3629
const X0 = [40 for i in 1:N_DAMS]
@@ -46,14 +39,9 @@ const B = [-1 0. 0. 0. 0. -1 0. 0. 0. 0.;
4639
0. 0. 1. -1 0. 0. 0. 1. -1 0.;
4740
0. 0. 0. 1. -1 0. 0. 0. 1. -1]
4841
# Define dynamic of the dam:
49-
function dynamic(t, x, u, w)
50-
return A*x + B*u + w
51-
end
52-
42+
dynamic(t, x, u, w) = A*x + B*u + w
5343
# Define cost corresponding to each timestep:
54-
function cost_t(t, x, u, w)
55-
return COST[t] * sum(u[1:N_DAMS])
56-
end
44+
cost_t(t, x, u, w) = COST[t] * sum(u[1:N_DAMS])
5745

5846
# We define here final cost a quadratic problem
5947
# we penalize the final costs if it is greater than 40.
@@ -65,18 +53,19 @@ function final_cost_dams(model, m)
6553
x = m[:x]
6654
u = m[:u]
6755
xf = m[:xf]
68-
@JuMP.variable(m, z1 >= 0)
69-
@JuMP.variable(m, z2 >= 0)
70-
@JuMP.variable(m, z3 >= 0)
71-
@JuMP.variable(m, z4 >= 0)
72-
@JuMP.variable(m, z5 >= 0)
73-
@JuMP.constraint(m, alpha == 0.)
74-
@JuMP.constraint(m, z1 >= 40 - xf[1])
75-
@JuMP.constraint(m, z2 >= 40 - xf[2])
76-
@JuMP.constraint(m, z3 >= 40 - xf[3])
77-
@JuMP.constraint(m, z4 >= 40 - xf[4])
78-
@JuMP.constraint(m, z5 >= 40 - xf[5])
79-
@JuMP.objective(m, Min, model.costFunctions(model.stageNumber-1, x, u, w) + 500. *(z1*z1+z2*z2+z3*z3+z4*z4+z5*z5))
56+
@JuMP.variable(m, z1 >= 0.0)
57+
@JuMP.variable(m, z2 >= 0.0)
58+
@JuMP.variable(m, z3 >= 0.0)
59+
@JuMP.variable(m, z4 >= 0.0)
60+
@JuMP.variable(m, z5 >= 0.0)
61+
@JuMP.constraint(m, alpha == 0.0)
62+
@JuMP.constraint(m, z1 >= 40.0 - xf[1])
63+
@JuMP.constraint(m, z2 >= 40.0 - xf[2])
64+
@JuMP.constraint(m, z3 >= 40.0 - xf[3])
65+
@JuMP.constraint(m, z4 >= 40.0 - xf[4])
66+
@JuMP.constraint(m, z5 >= 40.0 - xf[5])
67+
@JuMP.objective(m, Min, model.costFunctions(model.stageNumber-1, x, u, w)
68+
+ 500.0 * (z1*z1+z2*z2+z3*z3+z4*z4+z5*z5))
8069
end
8170

8271
##################################################
@@ -86,7 +75,7 @@ end
8675
const FORWARD_PASS = 10.
8776
const EPSILON = .01
8877
# Maximum number of iterations
89-
const MAX_ITER = 40
78+
const MAX_ITER = 10
9079
##################################################
9180

9281
"""Build probability distribution at each timestep.
@@ -113,13 +102,10 @@ function init_problem()
113102
X0, cost_t,
114103
dynamic, aleas,
115104
Vfinal=final_cost_dams)
116-
117105
# Add bounds for stocks:
118106
set_state_bounds(model, x_bounds)
119-
120-
121107
params = SDDPparameters(OPTIMIZER,
122-
passnumber=FORWARD_PASS,
108+
passnumber=1,
123109
compute_ub=10,
124110
gap=EPSILON,
125111
max_iterations=MAX_ITER)

‎examples/solver.jl‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,11 @@
11

2-
3-
#= using Clp =#
4-
#= SOLVER = ClpSolver() =#
5-
using Gurobi
2+
using Clp
3+
OPTIMIZER = optimizer_with_attributes(Clp.Optimizer, "LogLevel"=>0)
4+
#= using Xpress =#
5+
#= OPTIMIZER = optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG"=>0) =#
6+
#= using Gurobi =#
67
#SOLVER = Gurobi.GurobiSolver(OutputFlag=false, Threads=1)
7-
OPTIMIZER = optimizer_with_attributes(Gurobi.Optimizer,
8-
"OutputFlag"=>0)
8+
#= OPTIMIZER = optimizer_with_attributes(Gurobi.Optimizer, =#
9+
#= "OutputFlag"=>0) =#
910
#= using CPLEX =#
1011
#= SOLVER = CPLEX.CplexSolver(CPX_PARAM_SIMDISPLAY=0, CPX_PARAM_BARDISPLAY=0) =#

‎src/SDDPoptimize.jl‎

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -238,7 +238,7 @@ function finalpass!(sddp::SDDPInterface)
238238
println("\n", "#"^60)
239239
println("SDDP CONVERGENCE")
240240
@printf("- Exact lower bound: %.4e [Gap < %.2f%s]\n",
241-
lwb, 100*(upb+tol-lwb)/lwb, '%')
241+
lwb, 100*(upb+tol-lwb)/(1.0 + abs(lwb)), '%')
242242
@printf("- Estimation of upper-bound: %.4e\n", upb)
243243
@printf("- Upper-bound's s.t.d: %.4e\n", σ)
244244
@printf("- Confidence interval (%d%s): [%.4e , %.4e]",
@@ -290,7 +290,6 @@ $(SIGNATURES)
290290
function build_terminal_cost!(model::SPModel, problem::JuMP.Model,
291291
Vt::PolyhedralFunction, verbosity::Int64=0)
292292
# if shape is PolyhedralFunction, build terminal cost with it:
293-
294293
alpha = problem[:alpha]
295294
xf = problem[:xf]
296295
t = model.stageNumber -1
@@ -348,13 +347,16 @@ function build_model(model, param, t,verbosity::Int64=0)
348347
nw = model.dimNoises
349348

350349
# define variables in JuMP:
351-
@variable(m, model.xlim[i,t][1] <= x[i=1:nx] <= model.xlim[i,t][2])
352-
@variable(m, model.xlim[i,t][1] <= xf[i=1:nx]<= model.xlim[i,t][2])
353-
@variable(m, model.ulim[i,t][1] <= u[i=1:nu] <= model.ulim[i,t][2])
350+
@variable(m, model.xlim[i, t][1] <= x[i=1:nx] <= model.xlim[i, t][2])
351+
@variable(m, model.xlim[i, t][1] <= xf[i=1:nx]<= model.xlim[i, t][2])
352+
@variable(m, model.ulim[i, t][1] <= u[i=1:nu] <= model.ulim[i, t][2])
354353
@variable(m, alpha)
355354

356355
@variable(m, w[1:nw] == 0)
357-
m.ext[:cons] = @constraint(m, state_constraint, x .== 0)
356+
# This workaround is far from optimal. We should replace this line
357+
# with ParameterJuMP.jl
358+
@variable(m, x_constant[1:nx])
359+
m.ext[:cons] = @constraint(m, state_constraint, x .== x_constant)
358360

359361
@constraint(m, xf .== model.dynamics(t, x, u, w))
360362

@@ -388,7 +390,7 @@ function build_model(model, param, t,verbosity::Int64=0)
388390
if model.IS_SMIP
389391
m.colCat[2*nx+1:2*nx+nu] = model.controlCat
390392
end
391-
(verbosity >5) && print(m)
393+
(verbosity > 5) && println(m)
392394
return m
393395
end
394396

@@ -411,7 +413,8 @@ function build_model_dh(model, param, t, verbosity::Int64=0)
411413
@variable(m, model.xlim[i,t][1] <= xf[i=1:nx, j=1:ns]<= model.xlim[i,t][2])
412414
@variable(m, alpha[1:ns])
413415

414-
m.ext[:cons] = @constraint(m, state_constraint, x .== 0)
416+
@variable(m, x_constant[1:nx])
417+
m.ext[:cons] = @constraint(m, state_constraint, x .== x_constant)
415418

416419
for j=1:ns
417420
@constraint(m, xf[:, j] .== model.dynamics(t, x, u, ξ[:, j]))

0 commit comments

Comments
 (0)