Update tests (#2)

* Fix test

* Add tests on iterative method

* Fast usage values

* Add test precalculated

* Add seed for random

* Change solver

* Change solver

* Update ci.yml

* Update ci.yml

* Update sddp.jl

* Change saving directory for cuts

* Save protos in dev dir

* Saving cuts in dev directory

* Fix ci

---------

Co-authored-by: Juliette-Gerbaux <[email protected]>

.github/workflows/ci.yml

@@ -24,6 +24,9 @@ jobs:
           python -m pip install --upgrade pip
           pip install -r requirements.txt -r requirements-dev.txt
 
+      - name: Install Julia dependencies
+        run: julia -e 'import Pkg; Pkg.add("SDDP");Pkg.add("Clp");Pkg.add("JuMP");Pkg.add("PythonCall")'  
+
       - name: Test
         run: |
           pytest
@@ -32,3 +35,5 @@ jobs:
         run: |
           python -m mypy src
 
+
+

src/multi_stock_bellman_value_calculation.py

@@ -31,6 +31,7 @@
 from tqdm import tqdm
 from pathlib import Path
 import juliacall
+import os as os
 
 jl = juliacall.Main
 
@@ -218,6 +219,7 @@ def initialize_antares_problems(
     verbose: bool = False,
     save_protos: bool = False,
     load_from_protos: bool = False,
+    saving_dir: Optional[str] = None,
 ) -> Dict[TimeScenarioIndex, AntaresProblem]:
     """
     Creates Instances of the Antares problem for every week / scenario
@@ -241,15 +243,19 @@ def initialize_antares_problems(
         week_range = tqdm(week_range, desc="Problem initialization", colour="Yellow")
     for week in week_range:
         for scenario in range(param.len_scenario):
-            proto_path = (
-                output_path
-                + f"/protos/problem-{param.name_scenario[scenario]}-{week+1}.pkl"
-            )
-            already_processed = Path(proto_path).is_file() and load_from_protos
+            if saving_dir is not None:
+                proto_path = (
+                    saving_dir
+                    + f"/problem-{param.name_scenario[scenario]}-{week+1}.pkl"
+                )
+                already_processed = Path(proto_path).is_file() and load_from_protos
+            else:
+                already_processed = False
             m = AntaresProblem(
                 scenario=scenario,
                 week=week,
                 path=output_path,
+                saving_directory=saving_dir,
                 itr=1,
                 name_solver=name_solver,
                 name_scenario=param.name_scenario[scenario],
@@ -507,6 +513,7 @@ def get_all_costs(
     list_models: Dict[TimeScenarioIndex, AntaresProblem],
     multi_stock_management: MultiStockManagement,
     controls_list: np.ndarray,
+    saving_dir: Optional[str] = None,
     verbose: bool = False,
     already_init: bool = False,
     keep_intermed_res: bool = False,
@@ -529,7 +536,9 @@ def get_all_costs(
     tot_iter = 0
     times = []
     n_reservoirs = len(multi_stock_management.dict_reservoirs)
-    filename = "get_all_costs_run.pkl"
+    if keep_intermed_res or already_init:
+        assert saving_dir is not None
+        filename = saving_dir + "/get_all_costs_run.pkl"
 
     # Initializing the n_weeks*n_scenarios*n_controls*n_stocks values to fill
     shape_controls = list(controls_list.shape)
@@ -582,12 +591,13 @@ def Lget_costs(
     output_path: str,
     name_solver: str,
     controls_list: np.ndarray,
+    saving_directory: str,
     verbose: bool = False,
     direct_bellman_calc: bool = True,
     load_from_protos: bool = False,
     prefix: str = "",
 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
-    filename = f"{prefix}get_all_costs_run_{output_path[-27:]}.pkl"
+    filename = f"{saving_directory}/{prefix}get_all_costs_run_{output_path.replace('/','_')[-27:]}.pkl"
 
     # Initializing the n_weeks*n_scenarios*n_controls(*n_stocks) values to fill
     shape_controls = list(controls_list.shape)
@@ -609,15 +619,18 @@ def Lget_costs(
     for week in week_range:
         if week >= week_start:
             for scenario in range(param.len_scenario):
+                if not (os.path.exists(saving_directory)):
+                    os.makedirs(saving_directory)
                 proto_path = (
-                    output_path
-                    + f"/protos/problem-{param.name_scenario[scenario]}-{week+1}.pkl"
+                    saving_directory
+                    + f"/problem-{param.name_scenario[scenario]}-{week+1}.pkl"
                 )
                 already_processed = load_from_protos and Path(proto_path).is_file()
                 m = AntaresProblem(
                     scenario=scenario,
                     week=week,
                     path=output_path,
+                    saving_directory=saving_directory,
                     itr=1,
                     name_solver=name_solver,
                     name_scenario=param.name_scenario[scenario],
@@ -978,7 +991,7 @@ def select_controls_to_explore(
     """
     # This only a first version:
     if not rng:
-        rng = rng = np.random.default_rng()
+        rng = rng = np.random.default_rng(seed=12345)
     n_reservoirs = len(multi_stock_management.dict_reservoirs)
     n_weeks, n_scenarios = param.len_week, param.len_scenario
     controls_to_explore = np.zeros((n_weeks, n_scenarios, 1, n_reservoirs))
@@ -1237,6 +1250,7 @@ def cutting_plane_method(
     nSteps_bellman: int,
     method: str,
     correlations: np.ndarray,
+    saving_dir: str,
     maxiter: Optional[int] = None,
     precision: float = 5e-2,
     interp_mode: bool = False,
@@ -1382,6 +1396,7 @@ def cutting_plane_method(
             param=param,
             multi_stock_management=multi_stock_management,
             controls_list=controls_list,
+            saving_directory=saving_dir,
             output_path=output_path,
             name_solver=name_solver,
             verbose=verbose,
@@ -1446,6 +1461,7 @@ def iter_bell_vals(
     starting_pt: np.ndarray,
     nSteps_bellman: int,
     method: str,
+    saving_dir: str,
     name_solver: str = "CLP",
     precision: float = 1e-2,
     maxiter: int = 2,
@@ -1534,6 +1550,7 @@ def iter_bell_vals(
         param=param,
         multi_stock_management=multi_stock_management,
         output_path=output_path,
+        saving_directory=saving_dir,
         name_solver=name_solver,
         controls_list=controls_list,
         load_from_protos=True,
@@ -1574,6 +1591,7 @@ def iter_bell_vals(
             name_solver=name_solver,
             starting_pt=starting_pt,
             costs_approx=costs_approx,
+            saving_dir=saving_dir,
             costs=costs,
             future_costs_approx=future_costs_approx,
             nSteps_bellman=nSteps_bellman,
@@ -1626,6 +1644,7 @@ def sddp_cutting_planes(
     costs_approx: LinearCostEstimator,
     costs: np.ndarray,
     level_init: np.ndarray,
+    saving_dir: str,
     normalization: Dict[str, float],
     maxiter: Optional[int] = None,
     precision: float = 1e-2,
@@ -1675,6 +1694,7 @@ def sddp_cutting_planes(
         param.len_scenario,
         julia_reservoirs,
         julia_capp,
+        saving_dir,
         normalization["euro"],
         normalization["energy"],
     )
@@ -1688,6 +1708,7 @@ def sddp_cutting_planes(
             param.len_scenario,
             julia_reservoirs,
             julia_capp,
+            saving_dir,
             normalization["euro"],
             normalization["energy"],
         )
@@ -1716,6 +1737,7 @@ def sddp_cutting_planes(
             multi_stock_management=multi_stock_management,
             controls_list=controls,
             output_path=output_path,
+            saving_directory=saving_dir,
             name_solver=name_solver,
             verbose=False,
             load_from_protos=True,
@@ -1763,6 +1785,7 @@ def iter_bell_vals_v2(
     output_path: str,
     n_controls_init: int,
     starting_pt: np.ndarray,
+    saving_dir: str,
     normalization: Dict[str, float],
     name_solver: str = "CLP",
     precision: float = 1e-2,
@@ -1800,6 +1823,7 @@ def iter_bell_vals_v2(
         param=param,
         multi_stock_management=multi_stock_management,
         output_path=output_path,
+        saving_directory=saving_dir,
         name_solver=name_solver,
         controls_list=controls_list,
         load_from_protos=False,
@@ -1816,6 +1840,7 @@ def iter_bell_vals_v2(
         output_path=output_path,
         name_solver=name_solver,
         costs_approx=costs_approx,
+        saving_dir=saving_dir,
         costs=costs,
         level_init=starting_pt,
         precision=precision,

src/optimization.py

@@ -52,6 +52,7 @@ def __init__(
         week: int,
         path: str,
         itr: int = 1,
+        saving_directory: Optional[str] = None,
         name_solver: str = "CLP",
         name_scenario: int = -1,
         already_processed: bool = False,
@@ -67,6 +68,8 @@ def __init__(
             Week considered
         path:str :
             Path where mps files are stored
+        path:str :
+            Path where intermediate result are stored
         itr:int :
             Antares iteration considered (Default value = 1)
 
@@ -85,18 +88,19 @@ def __init__(
             mps_path = path + f"/problem-{name_scenario}-{week+1}--optim-nb-{itr}.mps"
             model = model_builder.ModelBuilder()  # type: ignore[no-untyped-call]
             model.import_from_mps_file(mps_path)
+            model_proto = model.export_to_proto()
+        else:
+            assert saving_directory is not None
+            proto_path = saving_directory + f"/problem-{name_scenario}-{week+1}.pkl"
             try:
-                model_proto = model.export_to_proto()
+                with open(proto_path, "rb") as file:
+                    model_proto, var_and_cstr_ids = pkl.load(file)
+                    self.stored_variables_and_constraints_ids = var_and_cstr_ids
             except FileNotFoundError:
                 print(
                     "Proto directory not found: Make sure the proto directory has been created within the mps directory"
                 )
                 raise FileNotFoundError
-        else:
-            proto_path = path + f"/protos/problem-{name_scenario}-{week+1}.pkl"
-            with open(proto_path, "rb") as file:
-                model_proto, var_and_cstr_ids = pkl.load(file)
-            self.stored_variables_and_constraints_ids = var_and_cstr_ids
 
         solver = pywraplp.Solver.CreateSolver(name_solver)
         assert solver, "Couldn't find any supported solver"

src/sddp.jl

@@ -1,11 +1,9 @@
 module Jl_SDDP
 
 using SDDP
-using HiGHS
 using Clp
 using JuMP
 using Base
-using Xpress
 using PythonCall
 
 # Define types for structured data
@@ -36,7 +34,7 @@ struct Normalizer
 end
 
 # Function to convert Python data to Julia structures
-function formater(n_weeks, n_scenarios, reservoirs_data, costs_approx_data, norm_euros::Float64=1e7, norm_enrgy::Float64=1e5)
+function formater(n_weeks, n_scenarios, reservoirs_data, costs_approx_data, saving_dir::String, norm_euros::Float64=1e7, norm_enrgy::Float64=1e5)
     round_energy = 4
     round_euro = 8
     round_price = 5
@@ -66,7 +64,7 @@ function formater(n_weeks, n_scenarios, reservoirs_data, costs_approx_data, norm
             round.(pyconvert(Matrix, costs_approx_data[i,j]["duals"])   /norm_price, digits=round_price), # € / MWh
             ) for j in 1:size_ca_data[2]] for i in 1:size_ca_data[1]]
 
-    return (n_weeks, n_scenarios, reservoirs, costs_approx, norms)
+    return (n_weeks, n_scenarios, reservoirs, costs_approx, norms, saving_dir)
 end
 
 
@@ -76,7 +74,7 @@ function generate_model(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Main.
         stages = 3*n_weeks,
         sense = :Min,
         lower_bound = 0.0,
-        optimizer = Xpress.Optimizer,
+        optimizer = Clp.Optimizer,
         # optimizer = HiGHS.Optimizer,
         # cut_oracle = SDDP.LevelOneCutOracle()
     ) do subproblem, stage
@@ -165,17 +163,17 @@ function stability_report(model)
     return SDDP.numerical_stability_report(model)
 end
 
-function reinit_cuts(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Main.Jl_SDDP.Reservoir}, costs_approx::Vector{Vector{Main.Jl_SDDP.LinInterp}}, norms::Normalizer)
+function reinit_cuts(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Main.Jl_SDDP.Reservoir}, costs_approx::Vector{Vector{Main.Jl_SDDP.LinInterp}}, norms::Normalizer, saving_dir::String)
     model = generate_model(n_weeks, n_scenarios, reservoirs, costs_approx, norms)
-    SDDP.write_cuts_to_file(model, "dev/cuts/sddp_current_cuts")
+    SDDP.write_cuts_to_file(model, saving_dir*"/sddp_current_cuts")
 end
 
-function manage_reservoirs(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Main.Jl_SDDP.Reservoir}, costs_approx::Vector{Vector{Main.Jl_SDDP.LinInterp}}, norms::Normalizer)
+function manage_reservoirs(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Main.Jl_SDDP.Reservoir}, costs_approx::Vector{Vector{Main.Jl_SDDP.LinInterp}}, norms::Normalizer, saving_dir::String)
     model = generate_model(n_weeks, n_scenarios, reservoirs, costs_approx, norms)
-    SDDP.read_cuts_from_file(model, "dev/cuts/sddp_current_cuts")
+    SDDP.read_cuts_from_file(model, saving_dir*"/sddp_current_cuts")
     # Training the model
     SDDP.train(model, stopping_rules = [SDDP.BoundStalling(40, 1e1)], iteration_limit = 700, cut_type = SDDP.MULTI_CUT)
-    SDDP.write_cuts_to_file(model, "dev/cuts/sddp_current_cuts")
+    SDDP.write_cuts_to_file(model, saving_dir*"/sddp_current_cuts")
 
     #Simulating
     simulation_results = get_trajectory(n_weeks, n_scenarios, reservoirs, model, norms)