Merge pull request #66 from leclere/integrate-ef

Integrate ef

Author: Henri-Gerard

examples/dynamic_example.jl

@@ -17,25 +17,21 @@ const SOLVER = ClpSolver()
 # const SOLVER = CplexSolver(CPX_PARAM_SIMDISPLAY=0)
 
 const N_STAGES = 5
-const N_SCENARIOS = 3
+const N_SCENARIOS = 1
 
-const DIM_STATES = 2
-const DIM_CONTROLS = 2
+const DIM_STATES = 1
+const DIM_CONTROLS = 1
 const DIM_ALEAS = 1
 
 
 
 #Constants that the user does not have to define
 const T0 = 1
 
-# Constants:
-const VOLUME_MAX = 100
-const VOLUME_MIN = 0
-
-const CONTROL_MAX = round(Int, .4/7. * VOLUME_MAX) + 1
+const CONTROL_MAX = round(Int, .4/7. * 100) + 1
 const CONTROL_MIN = 0
 
-const W_MAX = round(Int, .5/7. * VOLUME_MAX)
+const W_MAX = round(Int, .5/7. * 100)
 const W_MIN = 0
 const DW = 1
 
@@ -46,9 +42,7 @@ const ALEAS = linspace(W_MIN, W_MAX, N_ALEAS)
 const EPSILON = .05
 const MAX_ITER = 20
 
-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])
-
-const X0 = [50, 50]
+const X0 = 50*ones(DIM_STATES)
 
 Ax=[]
 Au=[]
@@ -146,23 +140,17 @@ function init_problem()
     x0 = X0
     aleas = generate_probability_laws()
 
-    x_bounds = [(VOLUME_MIN, VOLUME_MAX), (VOLUME_MIN, VOLUME_MAX)]
-    u_bounds  = []
-    for u = 1:DIM_CONTROLS
-        u_bounds = push!(u_bounds, (CONTROL_MIN, CONTROL_MAX))
-    end
- 
+    #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)
 
-    #set_state_bounds(model, x_bounds)
-
-    solver = SOLVER
-    params = SDDPparameters(solver, N_SCENARIOS, EPSILON, MAX_ITER)
+    params = SDDPparameters(SOLVER, N_SCENARIOS, EPSILON, MAX_ITER)
 
     return model, params
 end
@@ -177,10 +165,7 @@ function solve_dams(model,params,display=false)
 
     V, pbs = solve_SDDP(model, params, display)
 
-    aleas = simulate_scenarios(model.noises,
-                              (model.stageNumber,
-                               params.forwardPassNumber,
-                               model.dimNoises))
+    aleas = simulate_scenarios(model.noises,params.forwardPassNumber)
 
     params.forwardPassNumber = 1
 
@@ -189,7 +174,7 @@ function solve_dams(model,params,display=false)
     return stocks, V
 end
 
-#Solve the problem and try nb_iter times to generate radom data in case of infeasibility
+#Solve the problem and try nb_iter times to generate random data in case of infeasibility
 unsolve = true
 sol = 0
 i = 0
@@ -219,5 +204,6 @@ if (unsolve)
 else
     a,b = solve_dams(modelbis,paramsbis)
     println("solution =",sol)
+    println("V0 = ", b[1].lambdas[1,:]*X0+b[1].betas[1])
 end
 

src/StochDynamicProgramming.jl

@@ -12,9 +12,10 @@ module StochDynamicProgramming
 
 using JuMP, Distributions
 
-export solve_SDDP, NoiseLaw, simulate_scenarios,
+export solve_SDDP,
+            NoiseLaw, simulate_scenarios,
         SDDPparameters, LinearDynamicLinearCostSPmodel, set_state_bounds,
-        PiecewiseLinearCostSPmodel,
+        PiecewiseLinearCostSPmodel, extensive_formulation,
         PolyhedralFunction, NextStep, forward_simulations,
         StochDynProgModel, SDPparameters, solve_DP,
         sdp_forward_simulation, sampling, get_control, get_bellman_value

src/extensiveFormulation.jl

@@ -3,7 +3,8 @@
 #  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/.
 #############################################################################
-#  Define the extensive formulation to check the results of small problems
+#  Define the extensive formulation to check the results of small problems.
+#  The problem is instantiate on a tree.
 #############################################################################
 
 """
@@ -12,17 +13,14 @@ Contruct the scenario tree and solve the problem with measurability constraints
 function extensive_formulation(model,
                                param)
 
-
     #Recover all the constant in the model or in param
     laws = model.noises
-    N_NOISES = laws[1].supportSize
 
     DIM_STATE = model.dimStates
     DIM_CONTROL = model.dimControls
 
     X_init = model.initialState
 
-
     T = model.stageNumber-1
 
     mod = Model(solver=param.solver)
@@ -41,22 +39,21 @@ function extensive_formulation(model,
     #At each node, we have as many variables as nodes
     @defVar(mod,  u[t=1:T,n=1:DIM_CONTROL*N[t+1]])
     @defVar(mod,  x[t=1:T+1,n=1:DIM_STATE*N[t]])
-    @defVar(mod,  c[t=1:T,n=1:N_NOISES*N[t]])
+    @defVar(mod,  c[t=1:T,n=1:laws[t].supportSize*N[t]])
 
 
     #Computes the total probability of each node from the conditional probabilities
-    proba = []
-    push!(proba, laws[1].proba)
+    proba    = Vector{typeof(laws[1].proba)}(T)
+    proba[1] = laws[1].proba
     for t = 2 : T
-        push!(proba, zeros(N[t+1]))
+        proba[t] = zeros(N[t+1])
         for j = 1 : N[t]
-            for k = 1 : N_NOISES
-                proba[t][N_NOISES*(j-1)+k] = laws[t].proba[k]*proba[t-1][j]
+            for k = 1 : laws[t].supportSize
+                proba[t][laws[t].supportSize*(j-1)+k] = laws[t].proba[k]*proba[t-1][j]
             end
         end
     end
 
-
     #Instantiate the problem creating dynamic constraint at each node
     for t = 1 : (T)
         for n = 1 : N[t]
@@ -91,7 +88,7 @@ function extensive_formulation(model,
 
 
     #Define the objective of the function
-    @setObjective(mod, Min, sum{ sum{proba[t][N_NOISES*(n-1)+k]*c[t,N_NOISES*(n-1)+k],k = 1:N_NOISES} , t = 1:T, n=1:div(N[t+1],N_NOISES)})
+    @setObjective(mod, Min, sum{ sum{proba[t][laws[t].supportSize*(n-1)+k]*c[t,laws[t].supportSize*(n-1)+k],k = 1:laws[t].supportSize} , t = 1:T, n=1:div(N[t+1],laws[t].supportSize)})
 
     status = solve(mod)