Merge pull request #60 from leclere/dev_speedups

Fixes + unit tests

Author: leclere

examples/benchmark.jl

@@ -328,7 +328,7 @@ function benchmark_sdp()
     timing = zeros(n_benchmark)
     for n in 1:n_benchmark
         tic()
-        V_sdp = solve_DP(modelSDP, paramsSDP,1);
+        V_sdp = solve_DP(modelSDP, paramsSDP,0);
         timing[n] = toq()
     end
     @show timing

examples/stock-example.jl

@@ -76,8 +76,8 @@ if run_sdp
     infoStruct = "HD" # noise at time t is known before taking the decision at time t
 
     paramSDP = SDPparameters(spmodel, stateSteps, controlSteps, infoStruct)
-    Vs = sdp_optimize(spmodel,paramSDP)
-    lb_sdp = StochDynamicProgramming.get_value(spmodel,paramSDP,Vs)
+    Vs = solve_DP(spmodel,paramSDP, 1)
+    lb_sdp = StochDynamicProgramming.get_bellman_value(spmodel,paramSDP,Vs)
     println("Value obtained by SDP: "*string(lb_sdp))
 end
 

src/SDPoptimize.jl

@@ -144,7 +144,7 @@ function build_sdpmodel_from_spmodel(model::SPModel)
     if isa(model,PiecewiseLinearCostSPmodel)||isa(model,LinearDynamicLinearCostSPmodel)
         function cons_fun(t,x,u,w)
             for i in 1:model.dimStates
-                if (x[i]<=model.xlim[i][1]) || (x[i]>=model.xlim[i][2])
+                if (x[i]<model.xlim[i][1]) || (x[i]>model.xlim[i][2])
                     return false
                 end
             end
@@ -518,13 +518,13 @@ function sdp_forward_simulation(model::SPModel,
     nb_scenarios = size(scenarios)[2]
 
     costs = zeros(nb_scenarios)
-    states = zeros(TF,nb_scenarios)
-    controls = zeros(TF-1,nb_scenarios)
+    states = zeros(TF,nb_scenarios,model.dimStates)
+    controls = zeros(TF-1,nb_scenarios,model.dimControls)
 
 
     for k = 1:nb_scenarios
-        costs[k], states[:,k], controls[:,k] = sdp_forward_single_simulation(SDPmodel,
-                  param,scenarios[:,k],model.initialState,value,display)
+        costs[k], states[:,k,:], controls[:,k,:] = sdp_forward_single_simulation(SDPmodel,
+                  param,scenarios[:,k],model.initialState,V,display)
     end
 
     return costs, states, controls

test/runtests.jl

@@ -432,6 +432,8 @@ facts("SDP algorithm") do
 
             @fact test_costs --> true
 
+            @fact convertedSDPmodel.constraints(1,x,u,w) --> true
+
         end
 
         context("Solve and simulate using SDP") do
@@ -445,6 +447,14 @@ facts("SDP algorithm") do
                                                                                                         aleas_scen, x0,
                                                                                                         V_sdp, true )
 
+
+            costs_sdp2, stocks_sdp2, controls_sdp2 = StochDynamicProgramming.sdp_forward_simulation(modelSDP,
+                                                                                                    paramsSDP,
+                                                                                                    aleas_scen,
+                                                                                                    V_sdp, true )
+
+            @fact costs_sdp2[1] --> costs_sdp
+
             x = x0
             V_sdp = solve_DP(modelSDP, paramsSDP, false);
             V_sdp2 = StochDynamicProgramming.sdp_solve_HD(modelSDP, paramsSDP, false);
@@ -461,6 +471,15 @@ facts("SDP algorithm") do
             @fact v1 --> v2
             @fact (v1<=v3) --> true
 
+            paramsSDP.infoStructure = "DH"
+            costs_sdp3, stocks_sdp3, controls_sdp3 = StochDynamicProgramming.sdp_forward_simulation(modelSDP,
+                                                                                                    paramsSDP,
+                                                                                                    aleas_scen,
+                                                                                                    V_sdp3, true )
+            paramsSDP.infoStructure = "HD"
+
+            @fact costs_sdp3[1]>=costs_sdp2[1] --> true
+
             a,b = StochDynamicProgramming.generate_grid(modelSDP, paramsSDP)
 
             x_bounds = modelSDP.xlim
@@ -484,11 +503,15 @@ facts("SDP algorithm") do
             state_ref = zeros(2)
             state_ref[1] = stocks_sdp[2,1,1]
             state_ref[2] = stocks_sdp[2,1,2]
+            w = [4]
 
             @fact_throws get_control(modelSDP,paramsSDP,V_sdp3, 1, x)
+            @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x, w)[1] >= CONTROL_MIN) --> true
+            @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x, w)[1] <= CONTROL_MAX) --> true
+
             paramsSDP.infoStructure = "DH"
             @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x)[1] >= CONTROL_MIN) --> true
-            @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x)[1] >= CONTROL_MIN) --> true
+            @fact (get_control(modelSDP,paramsSDP,V_sdp3, 1, x)[1] <= CONTROL_MAX) --> true
 
             @fact size(stocks_sdp) --> (3,1,2)
             @fact size(controls_sdp) --> (2,1,2)