Merge pull request #111 from JuliaOpt/fix_sdp_bad_constraint_handling

Fix sdp bad constraint handling

Author: trigaut

examples/battery_storage_parallel.jl

@@ -77,7 +77,7 @@ println("library loaded")
     end
 
     function constraint(t, x, u, xi)
-    	return( (x[1] <= s_bounds[1][2] )&(x[1] >= s_bounds[1][1]))
+    	return true
     end
 
     function finalCostFunction(x)

examples/benchmark.jl

@@ -259,7 +259,7 @@ function benchmark_sdp()
     end
 
     function constraints(t, x, u, w)
-        return (VOLUME_MIN<=x[1]<=VOLUME_MAX)&(VOLUME_MIN<=x[2]<=VOLUME_MAX)
+        return true
     end
 
     function finalCostFunction(x)

src/SDPoptimize.jl

@@ -78,11 +78,6 @@ 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])
-                    return false
-                end
-            end
             return true
         end
         if in(:finalCostFunction,fieldnames(model))
@@ -418,7 +413,7 @@ function get_control(model::SPModel,param::SDPparameters,V, t::Int64, x::Array)
 
             next_state = SDPmodel.dynamics(t, x, u, w_sample)
 
-            if SDPmodel.constraints(t, next_state, u, w_sample)
+            if SDPmodel.constraints(t, x, u, w_sample)&&SDPutils.is_next_state_feasible(next_state, model.dimStates, model.xlim)
                 ind_next_state = SDPutils.real_index_from_variable(next_state, x_bounds, x_steps)
                 next_V = Vitp[ind_next_state...]
                 current_V += proba *(SDPmodel.costFunctions(t, x, u, w_sample) + next_V)
@@ -481,7 +476,7 @@ function get_control(model::SPModel,param::SDPparameters,V, t::Int64, x::Array,
 
         next_state = SDPmodel.dynamics(t, x, u, w)
 
-        if SDPmodel.constraints(t, next_state, u, w)
+        if SDPmodel.constraints(t, x, u, w)&&SDPutils.is_next_state_feasible(next_state, model.dimStates, model.xlim)
             ind_next_state = SDPutils.real_index_from_variable(next_state, x_bounds, x_steps)
             next_V = Vitp[ind_next_state...]
             current_V = SDPmodel.costFunctions(t, x, u, w) + next_V
@@ -586,7 +581,7 @@ function sdp_forward_single_simulation(model::StochDynProgModel,
 
                     next_state = model.dynamics(t, x, u, w_sample)
 
-                    if model.constraints(t, next_state, u, w_sample)
+                    if model.constraints(t, x, u, w_sample)&&SDPutils.is_next_state_feasible(next_state, model.dimStates, model.xlim)
                         ind_next_state = SDPutils.real_index_from_variable(next_state, x_bounds, x_steps)
                         next_V = Vitp[ind_next_state...]
                         current_V += proba *(model.costFunctions(t, x, u, w_sample) + next_V)
@@ -629,7 +624,7 @@ function sdp_forward_single_simulation(model::StochDynProgModel,
 
                 next_state = model.dynamics(t, x, u, scenario[t,1,:])
 
-                if model.constraints(t, next_state, u, scenario[t])
+                if model.constraints(t, x, u, scenario[t])&&SDPutils.is_next_state_feasible(next_state, model.dimStates, model.xlim)
                     ind_next_state = SDPutils.real_index_from_variable(next_state, x_bounds, x_steps)
                     next_V = Vitp[ind_next_state...]
                     current_V = model.costFunctions(t, x, u, scenario[t,1,:]) + next_V

src/SDPutils.jl

@@ -49,6 +49,34 @@ function real_index_from_variable(variable, bounds::Array, variable_steps::Array
     return tuple([1 + ( variable[i] - bounds[i][1] )/variable_steps[i] for i in 1:length(variable)]...)
 end
 
+"""
+Check if next state x_{t+1} satisfies state bounds constraints
+
+# Parameters
+* `next_stae::Array`:
+    the state we want to check
+* `x_dim::Int`:
+    the number of state variables
+* `x_bounds::Array`:
+    the state variables bounds
+
+# Returns
+* `index::Tuple{Float64}`:
+    the indexes of the variable
+
+"""
+function is_next_state_feasible(next_state, x_dim, x_bounds)
+
+    next_state_box_const = true
+
+    for i in 1:x_dim
+        next_state_box_const =  (next_state_box_const&&
+                                (next_state[i]>=x_bounds[i][1])&&
+                                (next_state[i]<=x_bounds[i][2]))
+    end
+
+    return next_state_box_const
+end
 
 """
 Computes the value function at time t evaluated at state x in a decision
@@ -104,7 +132,7 @@ function compute_V_given_x_t_DH(sampling_size, samples, probas, u_bounds,
             proba = probas[w]
             next_state = dynamics(t, x, u, w_sample)
 
-            if constraints(t, next_state, u, w_sample)
+            if constraints(t, x, u, w_sample)&&is_next_state_feasible(next_state, x_dim, x_bounds)
 
                 count_admissible_w = count_admissible_w + proba
                 ind_next_state = real_index_from_variable(next_state, x_bounds,
@@ -195,7 +223,8 @@ function compute_V_given_x_t_HD(sampling_size, samples, probas, u_bounds,
 
             next_state = dynamics(t, x, u, w_sample)
 
-            if constraints(t, next_state, u, w_sample)
+
+            if constraints(t, x, u, w_sample)&&is_next_state_feasible(next_state, x_dim, x_bounds)
                 admissible_u_w_count += 1
                 current_cost = cost(t, x, u, w_sample)
                 ind_next_state = real_index_from_variable(next_state, x_bounds,

test/sdp.jl

@@ -14,12 +14,16 @@ facts("Indexation for SDP") do
     ind = SDPutils.index_from_variable(var, bounds, steps)
     ind2 = SDPutils.real_index_from_variable(vart, bounds, steps)
 
+    checkFalse = SDPutils.is_next_state_feasible([0,1,2],3,bounds)
+    checkTrue = SDPutils.is_next_state_feasible([0.12,1.3,1.3],3,bounds)
+
 
     @fact ind --> (4,51,141)
     @fact ind2[1] --> roughly(4.2)
     @fact ind2[2] --> roughly(52.6)
     @fact ind2[3] --> roughly(144.2)
-
+    @fact checkFalse --> false
+    @fact checkTrue --> true
 
 end
 
@@ -67,7 +71,7 @@ facts("SDP algorithm") do
     end
 
     function constraints(t, x, u, w)
-        return (x[1]<=VOLUME_MAX)&(x[1]>=VOLUME_MIN)&(x[2]<=VOLUME_MAX)&(x[2]>=VOLUME_MIN)
+        return true
     end
 
     function finalCostFunction(x)