Merge pull request #115 from JuliaOpt/smip-sddp

Stochastic Mixed-Integer SDDP

Author: leclere

examples/benchmark.jl

@@ -167,20 +167,23 @@ function init_problem()
     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 = LinearDynamicLinearCostSPmodel(N_STAGES,
-                                                u_bounds,
-                                                x0,
-                                                cost_t,
-                                                dynamic,
-                                                aleas)
+    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, N_SCENARIOS, EPSILON, MAX_ITER)
+    params = SDDPparameters(solver,
+                            passnumber=N_SCENARIOS,
+                            gap=EPSILON,
+                            max_iterations=MAX_ITER)
 
     return model, params
 end

examples/dam.jl

@@ -142,16 +142,19 @@ function init_problem()
 
     x_bounds = [(0, 100)]
     u_bounds = [(0, 7), (0, 7)]
-    model = StochDynamicProgramming.LinearDynamicLinearCostSPmodel(N_STAGES,
-                                                u_bounds,
-                                                x0,
-                                                cost_t,
-                                                dynamic, aleas)
+    model = StochDynamicProgramming.LinearSPModel(N_STAGES,
+                                                  u_bounds,
+                                                  x0,
+                                                  cost_t,
+                                                  dynamic, aleas)
 
     set_state_bounds(model, x_bounds)
 
     solver = SOLVER
-    params = StochDynamicProgramming.SDDPparameters(solver, N_SCENARIOS, EPSILON, MAX_ITER)
+    params = StochDynamicProgramming.SDDPparameters(solver,
+                                                    passnumber=N_SCENARIOS,
+                                                    gap=EPSILON,
+                                                    max_iterations=MAX_ITER)
 
     return model, params
 end

examples/damsvalley.jl

@@ -107,21 +107,21 @@ function init_problem()
 
     x_bounds = [(VOLUME_MIN, VOLUME_MAX) for i in 1:N_DAMS]
     u_bounds = vcat([(CONTROL_MIN, CONTROL_MAX) for i in 1:N_DAMS], [(0., 200) for i in 1:N_DAMS]);
-    model = LinearDynamicLinearCostSPmodel(N_STAGES,
-                                                u_bounds,
-                                                X0,
-                                                cost_t,
-                                                dynamic,
-                                                aleas,
-                                                final_cost_dams)
+    model = LinearSPModel(N_STAGES, u_bounds,
+                          X0, cost_t,
+                          dynamic, aleas,
+                          final_cost_dams)
 
     # 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, FORWARD_PASS, EPSILON, MAX_ITER)
 
+    params = SDDPparameters(solver,
+                            passnumber=FORWARD_PASS,
+                            gap=EPSILON,
+                            max_iterations=MAX_ITER)
     return model, params
 end
 

examples/stock-example.jl

@@ -58,14 +58,18 @@ end
 ######## Setting up the SPmodel
 s_bounds = [(0, 1)] # bounds on the state
 u_bounds = [(CONTROL_MIN, CONTROL_MAX)] # bounds on controls
-spmodel = LinearDynamicLinearCostSPmodel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws)
+spmodel = LinearSPModel(N_STAGES,u_bounds,[S0],cost_t,dynamic,xi_laws)
 set_state_bounds(spmodel, s_bounds) # adding the bounds to the model
 println("Model set up")
 
 ######### Solving the problem via SDDP
 if run_sddp
     println("Starting resolution by SDDP")
-    paramSDDP = SDDPparameters(SOLVER, 10, 0, MAX_ITER) # 10 forward pass, stop at MAX_ITER
+    # 10 forward pass, stop at MAX_ITER
+    paramSDDP = StochDynamicProgramming.SDDPparameters(SOLVER,
+                                                    passnumber=10,
+                                                    gap=0,
+                                                    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)
     println("Lower bound obtained by SDDP: "*string(round(lb_sddp,4)))

src/SDDPoptimize.jl

@@ -37,6 +37,9 @@ fulfilled.
 
 """
 function solve_SDDP(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
     # initialize value functions:
     V, problems = initialize_value_functions(model, param)
     (verbose > 0) && println("Initial value function loaded into memory.")
@@ -47,6 +50,9 @@ end
 
 
 function solve_SDDP(model::SPModel, param::SDDPparameters, V::Vector{PolyhedralFunction}, verbose=0::Int64)
+    if model.IS_SMIP && isa(param.MIPSOLVER, Void)
+        error("MIP solver is not defined. Please set `param.MIPSOLVER`")
+    end
     # First step: process value functions if hotstart is called
     problems = hotstart_SDDP(model, param, V)
     sddp_stats = run_SDDP!(model, param, V, problems, verbose)
@@ -191,11 +197,8 @@ function estimate_upper_bound(model::SPModel, param::SDDPparameters,
                                 problem::Vector{JuMP.Model},
                                 n_simulation=1000::Int)
     aleas = simulate_scenarios(model.noises, n_simulation)
-    costs, stockTrajectories, _ = forward_simulations(model, param, problem, aleas)
-    return upper_bound(costs), costs
+    return estimate_upper_bound(model, param, aleas, problem)
 end
-
-
 function estimate_upper_bound(model::SPModel, param::SDDPparameters,
                                 aleas::Array{Float64, 3},
                                 problem::Vector{JuMP.Model})
@@ -204,12 +207,6 @@ function estimate_upper_bound(model::SPModel, param::SDDPparameters,
 end
 
 
-"""Build a collection of cuts initialized at 0"""
-function get_null_value_functions_array(model::SPModel)
-    return [PolyhedralFunction(zeros(1), zeros(1, model.dimStates), 1) for i in 1:model.stageNumber]
-end
-
-
 """
 Build a cut approximating terminal cost with null function
 
@@ -257,33 +254,35 @@ function build_models(model::SPModel, param::SDDPparameters)
 end
 
 function build_model(model, param, t)
-    m = Model(solver=param.solver)
+    m = Model(solver=param.SOLVER)
 
     nx = model.dimStates
     nu = model.dimControls
     nw = model.dimNoises
 
+    # define variables in JuMP:
     @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]<= model.xlim[i][2])
+    @variable(m,  model.ulim[i][1] <= u[i=1:nu] <=  model.ulim[i][2])
     @variable(m, alpha)
 
     @variable(m, w[1:nw] == 0)
     m.ext[:cons] = @constraint(m, state_constraint, x .== 0)
 
     @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)
     end
     if model.inequalityConstraints != nothing
         @constraint(m, model.inequalityConstraints(t, x, u, w) .<= 0)
     end
 
-    if typeof(model) == LinearDynamicLinearCostSPmodel
+    # Define objective function (could be linear or piecewise linear)
+    if isa(model.costFunctions, Function)
         @objective(m, Min, model.costFunctions(t, x, u, w) + alpha)
-
-    elseif typeof(model) == PiecewiseLinearCostSPmodel
+    elseif isa(model.costFunctions, Vector{Function})
         @variable(m, cost)
 
         for i in 1:length(model.costFunctions)
@@ -292,6 +291,11 @@ function build_model(model, param, t)
         @objective(m, Min, cost + alpha)
     end
 
+    # Add binary variable if problem is a SMIP:
+    if model.IS_SMIP
+        m.colCat[2*nx+1:2*nx+nu] = model.controlCat
+    end
+
     return m
 end
 
@@ -319,7 +323,7 @@ function initialize_value_functions(model::SPModel,
 
     solverProblems = build_models(model, param)
     V = PolyhedralFunction[
-                PolyhedralFunction([], Array{Float64}(0, model.dimStates), 0) for i in 1:model.stageNumber]
+                PolyhedralFunction(model.dimStates) for i in 1:model.stageNumber]
 
     # Build scenarios according to distribution laws:
     aleas = simulate_scenarios(model.noises, param.forwardPassNumber)
@@ -399,7 +403,7 @@ Bellman value (Float64)
 function get_bellman_value(model::SPModel, param::SDDPparameters,
                            t::Int64, Vt::PolyhedralFunction, xt::Vector{Float64})
 
-    m = Model(solver=param.solver)
+    m = Model(solver=param.SOLVER)
     @variable(m, alpha)
 
     for i in 1:Vt.numCuts
@@ -515,7 +519,7 @@ function exact_prune_cuts(model::SPModel, params::SDDPparameters, V::PolyhedralF
     ncuts = V.numCuts
     # Find all active cuts:
     if ncuts > 1
-        active_cuts = Bool[is_cut_relevant(model, i, V, params.solver) for i=1:ncuts]
+        active_cuts = Bool[is_cut_relevant(model, i, V, params.SOLVER) for i=1:ncuts]
         return PolyhedralFunction(V.betas[active_cuts], V.lambdas[active_cuts, :], sum(active_cuts))
     else
         return V

src/SDPoptimize.jl

@@ -76,7 +76,7 @@ function build_sdpmodel_from_spmodel(model::SPModel)
         return 0
     end
 
-    if isa(model,PiecewiseLinearCostSPmodel)||isa(model,LinearDynamicLinearCostSPmodel)
+    if isa(model,LinearSPModel)
         function cons_fun(t,x,u,w)
             return true
         end

src/StochDynamicProgramming.jl

@@ -10,12 +10,12 @@
 module StochDynamicProgramming
 include("SDPutils.jl")
 
-using JuMP, Distributions
+using MathProgBase, JuMP, Distributions
 
 export solve_SDDP,
         NoiseLaw, simulate_scenarios,
-        SDDPparameters, LinearDynamicLinearCostSPmodel, set_state_bounds,
-        PiecewiseLinearCostSPmodel, extensive_formulation,
+        SDDPparameters, LinearSPModel, set_state_bounds,
+         extensive_formulation,
         PolyhedralFunction, NextStep, forward_simulations,
         StochDynProgModel, SDPparameters, solve_DP,
         sdp_forward_simulation, sampling, get_control, get_bellman_value,

src/extensiveFormulation.jl

@@ -29,7 +29,7 @@ function extensive_formulation(model, param; verbose=0)
 
     X_init = model.initialState
     T = model.stageNumber-1
-    mod = Model(solver=param.solver)
+    mod = Model(solver=param.SOLVER)
 
     #Calculate the number of nodes n at each step on the scenario tree
     N = Array{Int64,2}(zeros(T+1,1))

src/forwardBackwardIterations.jl

@@ -213,10 +213,13 @@ function backward_pass!(model::SPModel,
                 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)
+                                                           t, state_t, alea_t,
+                                                           relaxation=model.IS_SMIP)
+
                 if solved
                     # We catch the subgradient λ:
                     subgradient_array[:, w] = nextstep.sub_gradient