Merge pull request #173 from mrchaos/master

Available with JuMP v0.19 or higher AND Julia 1.3.1 or higher available.

Author: frapac

.travis.yml

@@ -2,13 +2,13 @@ language: julia
 os:
   - linux
 julia:
-  - 0.6
+  - 1.3
+  - 1.4
 notifications:
   email: false
 
-script:
- - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
- - julia -e 'Pkg.clone(pwd()); Pkg.build("StochDynamicProgramming"); Pkg.test("StochDynamicProgramming"; coverage=true)'
 after_success:
-  - echo $TRAVIS_JULIA_VERSION
-  - julia -e 'Pkg.add("Coverage"); cd(Pkg.dir("StochDynamicProgramming")); using Coverage; Coveralls.submit(Coveralls.process_folder()); Codecov.submit(Codecov.process_folder())'
+  # push coverage results to Coveralls
+  - julia -e 'using Pkg; cd(Pkg.dir("StochDynamicProgramming")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
+  # push coverage results to Codecov
+  - julia -e 'using Pkg; cd(Pkg.dir("StochDynamicProgramming")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'

Project.toml

@@ -0,0 +1,32 @@
+name = "StochDynamicProgramming"
+uuid = "f6024387-75e5-5fdf-b095-327b05c0905f"
+authors = ["mrchaos <[email protected]>"]
+version = "0.5.0"
+
+[deps]
+CutPruners = "65d46eb8-70e9-5a30-bf48-2afa3a021b8f"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Nullables = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
+Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
+
+[compat]
+julia = "^1.3.0"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]

REQUIRE

@@ -31,7 +31,10 @@
 
 
 
-import StochDynamicProgramming, Distributions
+using StochDynamicProgramming, Distributions
+using Statistics
+using Distributed
+
 println("library loaded")
 
 # We have to define the instance on all the workers (processes)
@@ -62,7 +65,7 @@ println("library loaded")
 
     # create law of noises
     proba = 1/N_XI*ones(N_XI) # uniform probabilities
-    xi_support = collect(linspace(XI_MIN,XI_MAX,N_XI))
+    xi_support = collect(range(XI_MIN,stop=XI_MAX,length=N_XI))
     xi_law = StochDynamicProgramming.NoiseLaw(xi_support, proba)
     xi_laws = StochDynamicProgramming.NoiseLaw[xi_law for t in 1:N_STAGES-1]
 
@@ -112,4 +115,3 @@ lb_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs)
 println("Value obtained by SDP: "*string(lb_sdp))
 costsdp, states, stocks = StochDynamicProgramming.forward_simulations(spmodel,paramSDP,Vs,scenarios)
 println(mean(costsdp))
-

examples/battery_storage_parallel.jl

@@ -22,7 +22,7 @@ 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
 
 
 # COST:
-const COST = -66*2.7*(1 + .5*(rand(N_STAGES) - .5))
+const COST = -66*2.7*(1 .+ .5*(rand(N_STAGES) .- .5))
 
 # Constants:
 const VOLUME_MAX = 100
@@ -42,7 +42,7 @@ const X0 = [90]
 
 # Define aleas' space:
 const N_ALEAS = Int(round(Int, (W_MAX - W_MIN) / DW + 1))
-const ALEAS = linspace(W_MIN, W_MAX, N_ALEAS)
+const ALEAS = range(W_MIN,stop=W_MAX,length=N_ALEAS)
 
 
 # Define dynamic of the dam:
@@ -60,7 +60,7 @@ end
 in advance."""
 function solve_determinist_problem()
     println(alea_year)
-    m = Model(solver=SOLVER)
+    m = JuMP.Model(OPTIMIZER)
 
     @variable(m,  0           <= x[1:N_STAGES]  <= 100)
     @variable(m,  0.          <= u[1:N_STAGES-1]  <= 7)
@@ -74,10 +74,10 @@ function solve_determinist_problem()
 
     @constraint(m, x[1] .==X0)
 
-    status = solve(m)
+    status = JuMP.optimize!(m)
     println(status)
-    println(getObjectiveValue(m))
-    return getValue(u), getValue(x)
+    println(JuMP.getobjectivevalue(m))
+    return  JuMP.value.(u), JuMP.value.(x)
 end
 
 
@@ -86,7 +86,7 @@ function build_aleas()
     aleas = zeros(N_ALEAS, N_STAGES)
 
     # take into account seasonality effects:
-    unorm_prob = linspace(1, N_ALEAS, N_ALEAS)
+    unorm_prob = range(1,stop=N_ALEAS,length=N_ALEAS)
     proba1 = unorm_prob / sum(unorm_prob)
     proba2 = proba1[N_ALEAS:-1:1]
 
@@ -149,8 +149,8 @@ function init_problem()
 
     set_state_bounds(model, x_bounds)
 
-    solver = SOLVER
-    params = StochDynamicProgramming.SDDPparameters(solver,
+    optimizer = OPTIMIZER
+    params = StochDynamicProgramming.SDDPparameters(optimizer,
                                                     passnumber=N_SCENARIOS,
                                                     gap=EPSILON,
                                                     max_iterations=MAX_ITER)

examples/dam.jl

@@ -4,9 +4,9 @@
 
 ##################################################
 # Set a seed for reproductability:
-srand(2713)
-
-using StochDynamicProgramming, JuMP
+using Random
+using StochDynamicProgramming, JuMP, LinearAlgebra
+Random.seed!(2713)
 
 include("solver.jl")
 ##################################################
@@ -23,7 +23,7 @@ const N_ALEAS = 10
 
 # Cost are negative as we sell the electricity produced by
 # dams (and we want to minimize our problem)
-const COST = -66*2.7*(1 + .5*(rand(N_STAGES) - .5))
+const COST = -66*2.7*(1 .+ .5*(rand(N_STAGES) .- .5))
 
 # Constants:
 const VOLUME_MAX = 80
@@ -36,7 +36,7 @@ const CONTROL_MIN = 0
 const X0 = [40 for i in 1:N_DAMS]
 
 # Dynamic of stocks:
-const A = eye(N_DAMS)
+const A = I(N_DAMS)
 # The problem has the following structure:
 # dam1 -> dam2 -> dam3 -> dam4 -> dam5
 # We need to define the corresponding dynamic:
@@ -76,7 +76,7 @@ function final_cost_dams(model, m)
     @JuMP.constraint(m, z3 >= 40 - xf[3])
     @JuMP.constraint(m, z4 >= 40 - xf[4])
     @JuMP.constraint(m, z5 >= 40 - xf[5])
-    @JuMP.objective(m, Min, model.costFunctions(model.stageNumber-1, x, u, w) + 500.*(z1*z1+z2*z2+z3*z3+z4*z4+z5*z5))
+    @JuMP.objective(m, Min, model.costFunctions(model.stageNumber-1, x, u, w) + 500. *(z1*z1+z2*z2+z3*z3+z4*z4+z5*z5))
 end
 
 ##################################################
@@ -92,13 +92,13 @@ const MAX_ITER = 40
 """Build probability distribution at each timestep.
 Return a Vector{NoiseLaw}"""
 function generate_probability_laws()
-    laws = Vector{NoiseLaw}(N_STAGES-1)
+    laws = Vector{NoiseLaw}([])
     # uniform probabilities:
     proba = 1/N_ALEAS*ones(N_ALEAS)
 
     for t=1:N_STAGES-1
         support = rand(0:9, N_DAMS, N_ALEAS)
-        laws[t] = NoiseLaw(support, proba)
+        push!(laws, NoiseLaw(support, proba))
     end
     return laws
 end
@@ -118,7 +118,7 @@ function init_problem()
     set_state_bounds(model, x_bounds)
 
 
-    params = SDDPparameters(SOLVER,
+    params = SDDPparameters(OPTIMIZER,
                             passnumber=FORWARD_PASS,
                             compute_ub=10,
                             gap=EPSILON,
@@ -129,4 +129,3 @@ end
 # Solve the problem:
 model, params = init_problem()
 sddp = @time solve_SDDP(model, params, 2, 1)
-

examples/damsvalley.jl

@@ -7,12 +7,14 @@
 # Source: Adrien Cassegrain
 #############################################################################
 
-#srand(2713)
+using Random
+#Random.seed!(2713)
 
 using StochDynamicProgramming, JuMP, Clp
 
 #Constant that the user have to define himself
-const SOLVER = ClpSolver()
+const OPTIMIZER = optimizer_with_attributes(Clp.Optimizer,
+                                            "LogLevel"=>0)
 # const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0)
 
 const N_STAGES = 3
@@ -36,7 +38,7 @@ const DW = 1
 
 # Define aleas' space:
 const N_ALEAS = Int(round(Int, (W_MAX - W_MIN) / DW + 1))
-const ALEAS = linspace(W_MIN, W_MAX, N_ALEAS)
+const ALEAS = range(W_MIN,stop=W_MAX,length=N_ALEAS)
 
 const EPSILON = .05
 const MAX_ITER = 20
@@ -86,7 +88,7 @@ function build_aleas()
     aleas = zeros(N_ALEAS, N_STAGES)
 
     # take into account seasonality effects:
-    unorm_prob = linspace(1, N_ALEAS, N_ALEAS)
+    unorm_prob = range(1,stop=N_ALEAS,length=N_ALEAS)
     proba1 = unorm_prob / sum(unorm_prob)
     proba2 = proba1[N_ALEAS:-1:1]
 
@@ -120,13 +122,13 @@ Return a Vector{NoiseLaw}"""
 function generate_probability_laws()
     aleas = build_scenarios(N_SCENARIOS, build_aleas())
 
-    laws = Vector{NoiseLaw}(N_STAGES)
+    laws = Vector{NoiseLaw}([])
 
     # uniform probabilities:
     proba = 1/N_SCENARIOS*ones(N_SCENARIOS)
 
     for t=1:N_STAGES
-        laws[t] = NoiseLaw(aleas[:, t], proba)
+        push!(laws, NoiseLaw(aleas[:, t], proba))
     end
 
     return laws
@@ -142,15 +144,17 @@ function init_problem()
     #Define bounds for the control
     u_bounds  = [(CONTROL_MIN, CONTROL_MAX) for i in 1:DIM_CONTROLS]
 
-    model = LinearDynamicLinearCostSPmodel(N_STAGES,
-                                                u_bounds,
-                                                x0,
-                                                cost_t,
-                                                dynamic,
-                                                aleas)
-
-    params = SDDPparameters(SOLVER, N_SCENARIOS, EPSILON, MAX_ITER)
-
+    model = LinearSPModel(N_STAGES,
+                            u_bounds,
+                            x0,
+                            cost_t,
+                            dynamic,
+                            aleas)
+
+    params = SDDPparameters(OPTIMIZER,
+    passnumber=N_SCENARIOS,
+    gap=EPSILON,
+    max_iterations=MAX_ITER)
     return model, params
 end
 
@@ -181,6 +185,7 @@ i = 0
 nb_iter = 10
 
 while i<nb_iter
+    global model,params,i
     sol, firstControl, status = extensive_formulation(model,params)
     if (status == :Optimal)
         unsolve = false
@@ -207,4 +212,3 @@ else
     println("firstControl =", firstControl)
     println("V0 = ", b[1].lambdas[1,:]*X0+b[1].betas[1])
 end
-

examples/dynamic_example.jl

@@ -12,7 +12,7 @@
 #         u^i_t choosen knowing xi_1 .. xi_t
 #############################################################################
 
-using StochDynamicProgramming, Clp
+using StochDynamicProgramming, Clp, JuMP
 println("library loaded")
 
 run_sddp = true # false if you don't want to run sddp
@@ -21,7 +21,8 @@ test_simulation = false # false if you don't want to test your strategies
 
 ######## Optimization parameters  ########
 # choose the LP solver used.
-const SOLVER = ClpSolver() 			   # require "using Clp"
+const OPTIMIZER = optimizer_with_attributes(Clp.Optimizer,
+                        "LogLevel"=>0) 			   # require "using Clp"
 #const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0) # require "using CPLEX"
 
 # convergence test
@@ -49,7 +50,7 @@ const S0 = [0.5 for i=1:N_STOCKS]     # initial stock
 
 # create law of noises
 proba = 1/N_XI*ones(N_XI) # uniform probabilities
-xi_support = collect(linspace(XI_MIN,XI_MAX,N_XI))
+xi_support = collect(range(XI_MIN,stop=XI_MAX,length=N_XI))
 xi_law = StochDynamicProgramming.noiselaw_product([NoiseLaw(xi_support, proba) for i=1:N_STOCKS]...)
 xi_laws = NoiseLaw[xi_law for t in 1:N_STAGES-1]
 
@@ -77,22 +78,23 @@ println("Model set up")
 
 ######### Solving the problem via SDDP
 if run_sddp
-    tic()
+    @time begin
     println("Starting resolution by SDDP")
     # 10 forward pass, stop at MAX_ITER
-    paramSDDP = SDDPparameters(SOLVER,
+    paramSDDP = SDDPparameters(OPTIMIZER,
                                passnumber=1,
                                max_iterations=MAX_ITER)
     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();
+    println("Lower bound obtained by SDDP: "*string(round(lb_sddp,digits=4)))
+    end
+    println();
 end
 
 ######### Solving the problem via Dynamic Programming
 if run_sdp
-    tic()
+    @time begin
     println("Starting resolution by SDP")
     stateSteps = [step for i=1:N_STOCKS] # discretization step of the state
     controlSteps = [step for i=1:N_STOCKS] # discretization step of the control
@@ -104,8 +106,9 @@ if run_sdp
 
     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();
+    println("Value obtained by SDP: "*string(round(value_sdp,digits=4)))
+    end
+    println();
 end
 
 ######### Comparing the solutions on simulated scenarios.
@@ -115,6 +118,5 @@ if run_sddp && run_sdp && test_simulation
     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))*"%")
+            *string(round(200*mean(costsdp-costsddp)/abs(mean(costsddp+costsdp)),digits=3))*"%")
 end
-

examples/multistock-example.jl

@@ -1,9 +1,9 @@
 
-import StochDynamicProgramming
-
+using StochDynamicProgramming
+using Distributed
 
 """
-Solve SDDP in parallel, dispatching both forward and backward passes to process, 
+Solve SDDP in parallel, dispatching both forward and backward passes to process,
 which is not the most standard parallelization of SDDP.
 
 # Arguments