Feature/outer loop migrate criterion computation 2 (#789)

follows #780 
- compute outer loop criterion in Benders inner loop
- new mechanism to store vars indices which are linked to patterns
(PositiveUnsuppliedEnergy::\***AreaName\***)
- Update Python drivers so the outer loop can be launched like other
Xpansion steps
- Yaml parser (yaml-cpp) 
...

---------

Co-authored-by: Thomas Bittar <[email protected]>
Co-authored-by: Jason Maréchal <[email protected]>

.gitignore

@@ -37,6 +37,7 @@ _build
 _install
 antares-deps
 include/google
+vcpkg_installed
 
 # Cmake generated files
 CMakeCache.txt

data_test/external_loop_test/lp/description.txt

@@ -0,0 +1,44 @@
+We describe a system with 4 nodes, on 1 timestep, 2 scenarios :
+
+
+    |------- N1
+    |
+    N0 ----- N2
+    |
+    |------- N3
+
+    On node 0, there is a generator with :
+        - Initial P_max = 1
+        - Prod cost = 1
+        - Investment cost = 20
+
+    On nodes N1, N2, N3, there is a demand that differs between scenarios :
+        - On N1 : Demand1 = [1, 3]
+        - On N2 : Demand2 = [3, 3]
+        - On N3 : Demand3 = [0, 3]
+
+    Over the whole system we have :
+        - Spillage cost = 1
+        - ENS cost = 10
+
+    There are also transmission costs on lines :
+        - On L01 : transmission cost = 2
+        - On L02 : transmission cost = 1
+        - On L03 : transmission cost = 3
+
+    As the investment cost is higher than the ENS cost, adding 1MW of capacity would cost 20 to reduce ENS cost only by 10, hence an increased overall cost. Therefore the optimal investment is not to add any capacity beyond the existing one. Given the transmission cost, all flow will get to node 2.
+
+    This use case is simple enough so that we can count the number of hours with ENS with respect to the investment. Then we could use this test to check the behavior of any heuristic that aims at reaching a given target of LOLE (not done here).
+
+    We add 1 hour of ENS if there is at least 0.1 MWh of unsupplied energy, the "optimal" situation leads to 5h of ENS hours overall, hence 2.5h in expectation.
+
+    Each invested capacity will first go to node 2, then to node 1 and finally to node 3 given the transmission costs.
+
+    Hence we can deduce the number expected hours of ENS given the total capacity as follows  (invested capacity is total capacity - 1 as P_max = 1 initially):
+        - 0 <= P_max <= 2.9 : 2.5h (1 in area 1, 1 in area 2, 0.5 in area 3)
+        - 2.9 <= P_max <= 3.9 : 1.5h (1 in area 1, 0 in area 2, 0.5 in area 3)
+        - 3.9 <= P_max <= 5.9 : 1h (0.5 in area 1, 0 in area 2, 0.5 in area 3)
+        - 5.9 <= P_max <= 8.9 : 0.5h (0 in area 1, 0 in area 2, 0.5 in area 3)
+        - 8.9 <= P_max : 0h (0 in area 1, 0 in area 2, 0 in area 3)
+
+    Now, to test the heuristic taking into account the reliability constraint we can easily deduce the expected result depending on the criterion. For the test case, we require that all areas have less than 0.5 hours with ENS, hence, the solution is to have a total capacity of 3.9 (i.e. invest 2.9 as P_max = 1 initially)

data_test/external_loop_test/lp/options.json

@@ -4,26 +4,26 @@
     "RELATIVE_GAP": 1e-06,
     "RELAXED_GAP": 1e-05,
     "AGGREGATION": false,
-    "OUTPUTROOT": "data_test/external_loop_test/lp/",
+    "OUTPUTROOT": "data_test/external_loop_test/lp",
     "TRACE": true,
     "SLAVE_WEIGHT": "CONSTANT",
-    "SLAVE_WEIGHT_VALUE": 3,
+    "SLAVE_WEIGHT_VALUE": 0.5,
     "MASTER_NAME": "master",
-    "STRUCTURE_FILE": "data_test/external_loop_test/lp/structure.txt",
-    "INPUTROOT": "data_test/external_loop_test/lp/",
+    "STRUCTURE_FILE": "structure.txt",
+    "INPUTROOT": "data_test/external_loop_test/lp",
     "CSV_NAME": "benders_output_trace",
     "BOUND_ALPHA": true,
     "SEPARATION_PARAM": 0.5,
     "BATCH_SIZE": 0,
-    "JSON_FILE":"data_test/external_loop_test/expansion/out.json",
-    "LAST_ITERATION_JSON_FILE":"data_test/external_loop_test/expansion/last_iteration.json",
+    "JSON_FILE": "data_test/external_loop_test/expansion/out.json",
+    "LAST_ITERATION_JSON_FILE": "data_test/external_loop_test/expansion/last_iteration.json",
     "MASTER_FORMULATION": "integer",
     "SOLVER_NAME": "XPRESS",
     "TIME_LIMIT": 1000000000000.0,
-    "LOG_LEVEL": 1,
+    "LOG_LEVEL": 0,
     "LAST_MASTER_MPS": "master_last_iteration",
     "LAST_MASTER_BASIS": "master_last_basis.bss",
-    "EXT_LOOP_CRITERION_VALUE": 3.0,
-    "EXT_LOOP_CRITERION_TOLERANCE": 1e-1,
-    "EXT_LOOP_CRITERION_COUNT_THRESHOLD": 1e-1
+    "DO_OUTER_LOOP": true,
+    "OUTER_LOOP_OPTION_FILE": "data_test/external_loop_test/lp/outer_loop_options.yml",
+    "OUTER_LOOP_NUMBER_OF_SCENARIOS": 2
 }

data_test/external_loop_test/lp/outer_loop_options.yml

@@ -0,0 +1,15 @@
+# critère d'arrêt de l'algo
+stopping_threshold: 1e-4
+# seuil
+criterion_count_threshold: 1e-1
+ # tolerance entre seuil et valeur calculée
+criterion_tolerance: 1e-5
+patterns:
+  - area: "N0"
+    criterion: 1
+  - area: "N1"
+    criterion: 1
+  - area: "N2"
+    criterion: 1
+  - area: "N3"
+    criterion: 1

src/cpp/benders/CMakeLists.txt

@@ -4,7 +4,7 @@ add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/benders_core")
 add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/logger")
 add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/benders_sequential")
 add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/benders_by_batch")
-add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/external_loop")
+add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/outer_loop")
 
 
 add_subdirectory ("${CMAKE_CURRENT_SOURCE_DIR}/benders_mpi")

src/cpp/benders/benders_by_batch/BendersByBatch.cpp

@@ -6,6 +6,7 @@
 #include <numeric>
 
 #include "BatchCollection.h"
+#include "CustomVector.h"
 #include "RandomBatchShuffler.h"
 #include "glog/logging.h"
 BendersByBatch::BendersByBatch(
@@ -46,6 +47,10 @@ void BendersByBatch::InitializeProblems() {
       problem_count++;
     }
   }
+
+  // if (Rank() == rank_0) {
+  SetSubproblemsVariablesIndex();
+  // }
   init_problems_ = false;
 }
 void BendersByBatch::BroadcastSingleSubpbCostsUnderApprox() {
@@ -61,7 +66,6 @@ void BendersByBatch::Run() {
   if (init_data_) {
     PreRunInitialization();
   } else {
-    // only ?
     _data.stop = false;
   }
 
@@ -151,6 +155,9 @@ void BendersByBatch::SeparationLoop() {
     SolveBatches();
 
     if (Rank() == rank_0) {
+      outer_loop_criterion_.push_back(_data.outer_loop_current_iteration_data.outer_loop_criterion);
+      // TODO
+      //  UpdateOuterLoopMaxCriterionArea();
       UpdateTrace();
       SaveCurrentBendersData();
     }
@@ -195,8 +202,10 @@ void BendersByBatch::SolveBatches() {
     const auto &batch_sub_problems = batch.sub_problem_names;
     double batch_subproblems_costs_contribution_in_gap_per_proc = 0;
     double batch_subproblems_costs_contribution_in_gap = 0;
+    std::vector<double> external_loop_criterion_current_batch = {};
     BuildCut(batch_sub_problems,
-             &batch_subproblems_costs_contribution_in_gap_per_proc);
+             &batch_subproblems_costs_contribution_in_gap_per_proc,
+             external_loop_criterion_current_batch);
     Reduce(batch_subproblems_costs_contribution_in_gap_per_proc,
            batch_subproblems_costs_contribution_in_gap, std::plus<double>(),
            rank_0);
@@ -206,6 +215,9 @@ void BendersByBatch::SolveBatches() {
       _data.number_of_subproblem_solved += batch_sub_problems.size();
       _data.cumulative_number_of_subproblem_solved += batch_sub_problems.size();
       remaining_epsilon_ -= batch_subproblems_costs_contribution_in_gap;
+      // TODO
+      // AddVectors<double>(_data.outer_loop_current_iteration_data.outer_loop_criterion,
+      //                    external_loop_criterion_current_batch);
     }
 
     BroadCast(remaining_epsilon_, rank_0);
@@ -222,7 +234,8 @@ void BendersByBatch::SolveBatches() {
  */
 void BendersByBatch::BuildCut(
     const std::vector<std::string> &batch_sub_problems,
-    double *batch_subproblems_costs_contribution_in_gap_per_proc) {
+    double *batch_subproblems_costs_contribution_in_gap_per_proc,
+    std::vector<double> &external_loop_criterion_current_batch) {
   SubProblemDataMap subproblem_data_map;
   Timer subproblems_timer_per_proc;
   GetSubproblemCut(subproblem_data_map, batch_sub_problems,
@@ -235,7 +248,10 @@ void BendersByBatch::BuildCut(
   misprice_ = global_misprice;
   Gather(subproblem_data_map, gathered_subproblem_map, rank_0);
   SetSubproblemsWalltime(subproblems_timer_per_proc.elapsed());
-
+  // if (Options().EXTERNAL_LOOP_OPTIONS.DO_OUTER_LOOP) {
+  //   external_loop_criterion_current_batch =
+  //       ComputeSubproblemsContributionToOuterLoopCriterion(subproblem_data_map);
+  // }
   for (const auto &subproblem_map : gathered_subproblem_map) {
     for (auto &&[sub_problem_name, subproblem_data] : subproblem_map) {
       SetSubproblemCost(GetSubproblemCost() + subproblem_data.subproblem_cost);
@@ -270,6 +286,7 @@ void BendersByBatch::GetSubproblemCut(
       worker->fix_to(_data.x_cut);
       worker->solve(subproblem_data.lpstatus, Options().OUTPUTROOT,
                     Options().LAST_MASTER_MPS + MPS_SUFFIX, _writer);
+      worker->get_solution(subproblem_data.solution);
       worker->get_value(subproblem_data.subproblem_cost);  // solution phi(x,s)
       worker->get_subgradient(
           subproblem_data.var_name_and_subgradient);  // dual pi_s

src/cpp/benders/benders_by_batch/include/BendersByBatch.h

@@ -13,8 +13,8 @@ class BendersByBatch : public BendersMpi {
                  std::shared_ptr<MathLoggerDriver> mathLoggerDriver);
   ~BendersByBatch() override = default;
   void Run() override;
-  void BuildCut(const std::vector<std::string> &batch_sub_problems,
-                double *sum);
+  void BuildCut(const std::vector<std::string> &batch_sub_problems, double *sum,
+                std::vector<double> &external_loop_criterion_current_batch);
   std::string BendersName() const override { return "Benders By Batch mpi"; }
 
  protected: