Repository: Grid2op/chronix2grid (language: python)

Issue:
distinct std_noise for each data_type
change on line 196, noise_gen_tils, to enable distinct std for each data-type (load, solar, wind...)

    # retrieve the forecast parameters
    if data_type=='temperature':
        hs_mins, hs, std_hs = get_forecast_parameters(forecasts_params, load_params) 
    else:
        hs_mins, hs, std_hs = get_forecast_parameters(forecasts_params, load_params, data_type)


Discussion / hints:

Fix patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
index 9867461..827a0cf 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -37,8 +37,6 @@ jobs:
       - run: sudo apt-get install -y coinor-cbc
       - run: pip install .
       - run: pip install -r requirements.txt
-      - run: pip install pypsa==0.17.0
-      - run: pip install numpy==1.18.3
       - run: pip freeze
       - run:
           name: Run basic tests
diff --git a/.gitignore b/.gitignore
index 03db4b7..106af84 100644
--- a/.gitignore
+++ b/.gitignore
@@ -177,4 +177,5 @@ tests/data/output/generation/case118_l2rpn_wcci/agent_results/
 !tests/data/output/generation/expected_case118_l2rpn_neurips_1x
 !tests/data/output/generation/expected_case118_l2rpn_neurips_1x_modifySlackBeforeChronixGeneration
 !tests/data/output/generation/expected_case118_l2rpn_neurips_1x_hydro
-getting_started/example/input/generation/case118_l2rpn_wcci_Scenario_0/
\ No newline at end of file
+getting_started/example/input/generation/case118_l2rpn_wcci_Scenario_0/
+test_fix_forecast.py
\ No newline at end of file
diff --git a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/PypsaEconomicDispatch.py b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/PypsaEconomicDispatch.py
index eb2b821..c4ed081 100644
--- a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/PypsaEconomicDispatch.py
+++ b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/PypsaEconomicDispatch.py
@@ -12,7 +12,14 @@
 import numpy as np
 import pandas as pd
 from collections import namedtuple
-import pypsa
+import warnings
+
+with warnings.catch_warnings():
+    warnings.filterwarnings("ignore", category=RuntimeWarning)
+    warnings.filterwarnings("ignore", category=DeprecationWarning)
+    # TODO logger !
+    # some warnings on some pypsa / numpy versions
+    import pypsa
 
 from ._EDispatch_L2RPN2020.run_economic_dispatch import main_run_disptach
 
@@ -120,7 +127,10 @@ def from_dataframe(cls, env_df):
         -------
         net: :class:`pypsa.Network`
         """
-        net = cls()
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", category=FutureWarning)
+            # issue with pypsa on some version of numpy 
+            net = cls()
         net._df = env_df
 
         carrier_types_to_exclude = ['wind', 'solar']
diff --git a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py
index 4b49ea0..c0f042d 100644
--- a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py
+++ b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py
@@ -9,7 +9,7 @@
 import numpy as np
 import pandas as pd
 import copy 
-import pypsa
+import warnings
 
 from chronix2grid.generation.dispatch.utils import RampMode
 
@@ -268,7 +268,7 @@ def get_grouped_snapshots(snapshot, mode):
     """    
     # Define all posibilities mode
     periods = {'day': snapshot.groupby(snapshot.day).values(),
-               'week': snapshot.groupby(snapshot.week).values(),
+               'week': snapshot.groupby(snapshot.isocalendar().week).values(),
                'month': snapshot.groupby(snapshot.month).values()
     }
     return periods[mode]
@@ -307,7 +307,7 @@ def run_opf(net,
         Results of OPF dispatch
     """    
     to_disp = {'day': demand.index.day.unique().values[0],
-               'week': demand.index.week.unique().values[0],
+               'week': demand.index.isocalendar().week.unique()[0],
                'month': demand.index.month.unique().values[0],
     }
     mode = params['mode_opf']
@@ -351,7 +351,10 @@ def run_opf(net,
     net.generators_t.p_min_pu = net.generators_t.p_min_pu.iloc[0:0, 0:0]
     
     # Set snapshots
-    net.set_snapshots(demand.index)
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", category=FutureWarning)
+        # pypsa and pandas version
+        net.set_snapshots(demand.index)
     
     
     # ++  ++  ++  ++  ++  ++  ++  ++  ++  ++  ++ 
diff --git a/chronix2grid/generation/consumption/ConsumptionGeneratorBackend.py b/chronix2grid/generation/consumption/ConsumptionGeneratorBackend.py
index 3462572..722ac1c 100644
--- a/chronix2grid/generation/consumption/ConsumptionGeneratorBackend.py
+++ b/chronix2grid/generation/consumption/ConsumptionGeneratorBackend.py
@@ -38,11 +38,11 @@ def __init__(self, out_path, seed, params, loads_charac, load_config_manager, wr
         self.out_path = out_path
         self.day_lag = day_lag
 
-    def run(self, load_weekly_pattern=None):
+    def run(self, load_weekly_pattern=None, return_ref_curve=False, use_legacy=True):
         """
         Runs the generation model in ``chronix2grid.generation.consumption.generate_load`` and writes chronics
         """
         if load_weekly_pattern is None:
             load_weekly_pattern = self.load_config_manager.read_specific()
         return main(self.out_path, self.seed, self.params, self.loads_charac, load_weekly_pattern, self.write_results,
-                    day_lag=self.day_lag)
+                    day_lag=self.day_lag, return_ref_curve=return_ref_curve, use_legacy=use_legacy)
diff --git a/chronix2grid/generation/consumption/consumption_utils.py b/chronix2grid/generation/consumption/consumption_utils.py
index 841e763..9f772e8 100755
--- a/chronix2grid/generation/consumption/consumption_utils.py
+++ b/chronix2grid/generation/consumption/consumption_utils.py
@@ -7,8 +7,6 @@
 # This file is part of Chronix2Grid, A python package to generate "en-masse" chronics for loads and productions (thermal, renewable)
 
 import os
-import calendar
-from datetime import datetime
 import numpy as np
 import pandas as pd
 from scipy.interpolate import interp1d
@@ -16,11 +14,21 @@
 from .. import generation_utils as utils
 import chronix2grid.constants as cst
 
-def compute_loads(loads_charac, temperature_noise, params, load_weekly_pattern, start_day, add_dim, day_lag=0):
+def compute_loads(loads_charac, temperature_noise, params, load_weekly_pattern,
+                  start_day, add_dim, day_lag=0,
+                  return_ref_curve=False,
+                  use_legacy=True):
     #6  # this is only TRUE if you simulate 2050 !!! formula does not really work
     
     # Compute active part of loads
     weekly_pattern = load_weekly_pattern['test'].values
+    if use_legacy:
+        isoweekday = None
+        hour_minutes = None
+    else:
+        datetime_lwp = pd.to_datetime(load_weekly_pattern["datetime"], format="%Y-%m-%d %H:%M:%S")
+        isoweekday = np.array([el.isoweekday() for el in datetime_lwp])
+        hour_minutes = np.array([el.hour * 60 + el.minute for el in datetime_lwp])
     
     # start_day_of_week = start_day.weekday()
     # first_dow_chronics = datetime.strptime(load_weekly_pattern["datetime"].iloc[1], "%Y-%m-%d %H:%M:%S").weekday()
@@ -28,20 +36,52 @@ def compute_loads(loads_charac, temperature_noise, params, load_weekly_pattern,
     # day_lag = (first_dow_chronics - start_day_of_week) % 7
     # day_lag = 0
     loads_series = {}
+    ref_curves = None
     for i, name in enumerate(loads_charac['name']):
         mask = (loads_charac['name'] == name)
         if loads_charac[mask]['type'].values == 'residential':
             locations = [loads_charac[mask]['x'].values[0], loads_charac[mask]['y'].values[0]]
             Pmax = loads_charac[mask]['Pmax'].values[0]
-            loads_series[name] = compute_residential(locations, Pmax, temperature_noise, params, weekly_pattern, index=i, day_lag=day_lag, add_dim=add_dim)
+            tmp_ = compute_residential(locations, Pmax, temperature_noise,
+                                       params, weekly_pattern, index=i,
+                                       day_lag=day_lag, add_dim=add_dim,
+                                       return_ref_curve=return_ref_curve,
+                                       isoweekday_lwp=isoweekday,
+                                       hour_minutes_lwp=hour_minutes)
+            if return_ref_curve:
+                if ref_curves is None:
+                    ref_curves = np.zeros((tmp_[0].shape[0], loads_charac.shape[0]))
+                loads_series[name], ref_curves[:,i] = tmp_
+            else:
+                loads_series[name] = tmp_
 
         if loads_charac[mask]['type'].values == 'industrial':
             raise NotImplementedError("Impossible to generate industrial loads for now.")
             Pmax = loads_charac[mask]['Pmax'].values[0]
             loads_series[name] = compute_industrial(Pmax, params)
-    return loads_series
+    if not return_ref_curve:
+        return loads_series
+    return loads_series, ref_curves
+
+
+def get_seasonal_pattern(params):
+    Nt_inter = int(params['T'] // params['dt'] + 1)
+    t = np.linspace(0., (params['end_date'] - params["start_date"]).total_seconds(), Nt_inter, endpoint=True, dtype=float)
+    start_year = pd.to_datetime(str(params['start_date'].year) + '-01-01', format='%Y-%m-%d')
+    start_min = float(pd.Timedelta(params['start_date'] - start_year).total_seconds())
+    nb_sec_per_day =  24. * 60. * 60.
+    nb_sec_per_year = (365. * nb_sec_per_day)
+    year_pattern = 2. * np.pi / nb_sec_per_year
+    seasonal_pattern = 1.5 / 7. * np.cos(year_pattern * (t + start_min - 45 * nb_sec_per_day))  # min of the load is 15 of February so 45 days after beginning of year
+    seasonal_pattern += 5.5 / 7.
+    return seasonal_pattern
+
 
-def compute_residential(locations, Pmax, temperature_noise, params, weekly_pattern, index, day_lag=None, add_dim=0):
+def compute_residential(locations, Pmax, temperature_noise, params,
+                        weekly_pattern, index, day_lag=None, add_dim=0,
+                        return_ref_curve=False,
+                        isoweekday_lwp=None,
+                        hour_minutes_lwp=None):
 
 
     # Compute refined signals
@@ -54,28 +94,18 @@ def compute_residential(locations, Pmax, temperature_noise, params, weekly_patte
     temperature_signal = temperature_signal.astype(float)
     
     # Compute seasonal pattern
-    Nt_inter = int(params['T'] // params['dt'] + 1)
-    
-    # t = np.linspace(0, params['T'], Nt_inter, endpoint=True)
-    t = np.linspace(0., (params['end_date'] - params["start_date"]).total_seconds(), Nt_inter, endpoint=True, dtype=float)
-    
-    start_year = pd.to_datetime(str(params['start_date'].year) + '/01/01', format='%Y-%m-%d')
-    start_min = float(pd.Timedelta(params['start_date'] - start_year).total_seconds())
-    nb_sec_per_day =  24. * 60. * 60.
-    nb_sec_per_year = (365. * nb_sec_per_day)
-    year_pattern = 2. * np.pi / nb_sec_per_year
-    seasonal_pattern = 1.5 / 7. * np.cos(year_pattern * (t + start_min - 45 * nb_sec_per_day))  # min of the load is 15 of February so 45 days after beginning of year
-    #seasonal_pattern = 1.5 / 7. * np.cos(year_pattern * (t - start_min - 30 * nb_sec_per_day)) #older version to be removed
-    seasonal_pattern += 5.5 / 7.
+    seasonal_pattern = get_seasonal_pattern(params)
 
     # Get weekly pattern
-    weekly_pattern = compute_load_pattern(params, weekly_pattern, index, day_lag)
+    weekly_pattern = compute_load_pattern(params, weekly_pattern, index, day_lag, isoweekday_lwp=isoweekday_lwp, hour_minutes_lwp=hour_minutes_lwp)
+    
     std_temperature_noise = params['std_temperature_noise']
     residential_series = Pmax * weekly_pattern * (std_temperature_noise * temperature_signal + seasonal_pattern)
-
+    if return_ref_curve:
+        return residential_series, Pmax * weekly_pattern * seasonal_pattern
     return residential_series
 
-def compute_load_pattern(params, weekly_pattern, index, day_lag):
+def compute_load_pattern(params, weekly_pattern, index, day_lag=None, isoweekday_lwp=None, hour_minutes_lwp=None):
     """
     Loads a typical hourly pattern, and interpolates it to generate
     a smooth solar generation pattern between 0 and 1
@@ -89,38 +119,57 @@ def compute_load_pattern(params, weekly_pattern, index, day_lag):
         (np.array) A smooth solar pattern
     """
     # solar_pattern resolution : 1H, 8761
-
-    if day_lag is None:
-        nb_step_lag_for_starting_day = 0
-    else:
-        nb_step_lag_for_starting_day = 12 * 24 * day_lag
-        
+    
     # Keep only one week of pattern
     index_weekly_perweek = 12 * 24 * 7
-    index %= int((weekly_pattern.shape[0] - nb_step_lag_for_starting_day) / index_weekly_perweek - 1)
-    weekly_pattern = weekly_pattern[(nb_step_lag_for_starting_day + index * index_weekly_perweek):(nb_step_lag_for_starting_day + (index + 1) * index_weekly_perweek)]
+    if isoweekday_lwp is None or hour_minutes_lwp is None:
+        # try to guess where to start from input data
+        if day_lag is None:
+            nb_step_lag_for_starting_day = 0
+        else:
+            nb_step_lag_for_starting_day = 12 * 24 * day_lag
+        index %= int((weekly_pattern.shape[0] - nb_step_lag_for_starting_day) / index_weekly_perweek - 1)
+        first_index = (nb_step_lag_for_starting_day + index * index_weekly_perweek)
+    else:
+        # be smarter and take a week starting the same weekday at the same hour than the params["start_date"]
+        isoweekday_start = params["start_date"].isoweekday()
+        iso_hm_start = params["start_date"].hour * 60 + params["start_date"].minute
+        possible_first_index = np.where((isoweekday_lwp == isoweekday_start) & (iso_hm_start == hour_minutes_lwp))[0]
+        index_modulo = index % possible_first_index.shape[0]
+        first_index = possible_first_index[index_modulo]
+        
+    # now extract right week of data
+    last_index = first_index + index_weekly_perweek
+    weekly_pattern = weekly_pattern[first_index:last_index]
     weekly_pattern /= np.mean(weekly_pattern)
 
-    start_year = pd.to_datetime(str(params['start_date'].year) + '/01/01', format='%Y-%m-%d')
-    T_bis = int(pd.Timedelta(params['end_date'] - start_year).total_seconds() // (60))
-
-    Nt_inter_hr = int(T_bis // 5 + 1)
-    N_repet = int((Nt_inter_hr - 1) // len(weekly_pattern) + 1)
-    stacked_weekly_pattern = weekly_pattern
-    for i in range(N_repet - 1):
-        stacked_weekly_pattern = np.append(stacked_weekly_pattern, weekly_pattern)
+    # now generate an ouput of the right length
+    if isoweekday_lwp is None or hour_minutes_lwp is None:
+        # legacy usage... does not work at all I don't know why
+        start_year = pd.to_datetime(str(params['start_date'].year) + '-01-01', format='%Y-%m-%d')
+        T_bis = int(pd.Timedelta(params['end_date'] - start_year).total_seconds() // (60))
 
-    # The time is in minutes
-    t_pattern = np.linspace(0, 60 * 7 * 24 * N_repet, 12 * 7 * 24 * N_repet, endpoint=False)
-    f2 = interp1d(t_pattern, stacked_weekly_pattern, kind='cubic')
-
-    Nt_inter = int(params['T'] // params['dt'] + 1)
-    start_year = pd.to_datetime(str(params['start_date'].year) + '/01/01', format='%Y-%m-%d')
-    start_min = int(pd.Timedelta(params['start_date'] - start_year).total_seconds() // 60)
-    end_min = int(pd.Timedelta(params['end_date'] - start_year).total_seconds() // 60)
-    t_inter = np.linspace(start_min, end_min, Nt_inter, endpoint=True)
-    output = f2(t_inter)
-    output = output * (output > 0)
+        Nt_inter_hr = int(T_bis // 5 + 1)
+        N_repet = int((Nt_inter_hr - 1) // len(weekly_pattern) + 1)
+        stacked_weekly_pattern = np.tile(weekly_pattern, N_repet)
+        
+        # The time is in minutes
+        t_pattern = np.linspace(0, 60 * 7 * 24 * N_repet, 12 * 7 * 24 * N_repet, endpoint=False)
+        f2 = interp1d(t_pattern, stacked_weekly_pattern, kind='cubic')
+
+        Nt_inter = int(params['T'] // params['dt'] + 1)
+        start_year = pd.to_datetime(str(params['start_date'].year) + '-01-01', format='%Y-%m-%d')
+        start_min = int(pd.Timedelta(params['start_date'] - start_year).total_seconds() // 60)
+        end_min = int(pd.Timedelta(params['end_date'] - start_year).total_seconds() // 60)
+        t_inter = np.linspace(start_min, end_min, Nt_inter, endpoint=True)
+        output = f2(t_inter)
+        output = output * (output > 0)
+    else:
+        # new usage
+        nb_ts = int((params['end_date'] - params['start_date']).total_seconds() / 60 / params["dt"] + 1)  # +1 is because of the buggy stuff above...
+        N_repet = np.ceil(nb_ts / weekly_pattern.shape[0]).astype(int)
+        stacked_weekly_pattern = np.tile(weekly_pattern, N_repet)
+        output = stacked_weekly_pattern[:nb_ts]
 
     return output
 
diff --git a/chronix2grid/generation/consumption/generate_load.py b/chronix2grid/generation/consumption/generate_load.py
index 483de97..e7ce0b1 100755
--- a/chronix2grid/generation/consumption/generate_load.py
+++ b/chronix2grid/generation/consumption/generate_load.py
@@ -19,7 +19,22 @@
 from .. import generation_utils as utils
 
 
-def main(scenario_destination_path, seed, params, loads_charac, load_weekly_pattern, write_results = True, day_lag=0):
+def get_add_dim(params, loads_charac):
+    add_dim = 0
+    dx_corr = int(params['dx_corr'])
+    dy_corr = int(params['dy_corr'])
+    for x,y  in zip(loads_charac["x"], loads_charac["y"]):
+        x_plus = int(x // dx_corr + 1)
+        y_plus = int(y // dy_corr + 1)
+        add_dim = max(y_plus, add_dim)
+        add_dim = max(x_plus, add_dim)
+    return add_dim
+
+
+def main(scenario_destination_path, seed, params, loads_charac,
+         load_weekly_pattern, write_results = True, day_lag=0,
+         return_ref_curve=False,
+         use_legacy=True):
     """
     This is the load generation function, it allows you to generate consumption chronics based on demand nodes characteristics and on weekly demand patterns.
 
@@ -48,16 +63,7 @@ def main(scenario_destination_path, seed, params, loads_charac, load_weekly_patt
         end=params['end_date'],
         freq=str(params['dt']) + 'min')
 
-
-    add_dim = 0
-    dx_corr = int(params['dx_corr'])
-    dy_corr = int(params['dy_corr'])
-    for x,y  in zip(loads_charac["x"], loads_charac["y"]):
-        x_plus = int(x // dx_corr + 1)
-        y_plus = int(y // dy_corr + 1)
-        add_dim = max(y_plus, add_dim)
-        add_dim = max(x_plus, add_dim)
-        #add_dim=0 #to get back to when this parameter did not exist - to be removed
+    add_dim = get_add_dim(params, loads_charac)
     
     # Generate GLOBAL temperature noise
     print('Computing global auto-correlated spatio-temporal noise for thermosensible demand...') ## temperature is simply to reflect the fact that loads is correlated spatially, and so is the real "temperature". It is not the real temperature.
@@ -71,7 +77,11 @@ def main(scenario_destination_path, seed, params, loads_charac, load_weekly_patt
                                        load_weekly_pattern,
                                        start_day=start_day,
                                        add_dim=add_dim,
-                                       day_lag=day_lag)
+                                       day_lag=day_lag,
+                                       return_ref_curve=return_ref_curve,
+                                       use_legacy=use_legacy)
+    if return_ref_curve:
+        loads_series, ref_curve = loads_series
     loads_series['datetime'] = datetime_index
 
     # Save files
@@ -91,5 +101,7 @@ def main(scenario_destination_path, seed, params, loads_charac, load_weekly_patt
         write_results=write_results,
         index=False
     )
-    
-    return load_p, load_p_forecasted
+    if not return_ref_curve:
+        return load_p, load_p_forecasted
+    else:
+        return load_p, load_p_forecasted, ref_curve
diff --git a/chronix2grid/generation/generation_utils.py b/chronix2grid/generation/generation_utils.py
index e92f4cb..1ddeef9 100644
--- a/chronix2grid/generation/generation_utils.py
+++ b/chronix2grid/generation/generation_utils.py
@@ -32,21 +32,7 @@ def make_generation_input_output_directories(input_folder, case, year, output_fo
     return dispatch_input_folder, dispatch_input_folder_case, dispatch_output_folder
 
 
-def generate_coarse_noise(prng, params, data_type, add_dim):
-    """
-    This function generates a spatially and temporally correlated noise.
-    Because it may take a lot of time to compute a correlated noise on
-    a too fine mesh, we recommend to first compute a correlated signal
-    on a coarse mesh, and then use the interpolation function.
-
-    Input:
-        params: (dict) Defines the mesh dimensions and
-            precision. Also define the correlation scales
-
-    Output:
-        (np.array) 3D autocorrelated noise
-    """
-
+def get_nx_ny_nt(data_type, params, add_dim):
     # Get computation domain size
     Lx = params['Lx']
     Ly = params['Ly']
@@ -61,7 +47,26 @@ def generate_coarse_noise(prng, params, data_type, add_dim):
     Nx_comp = int(Lx // dx_corr + 1) + add_dim
     Ny_comp = int(Ly // dy_corr + 1) + add_dim
     Nt_comp = int(T // dt_corr + 1) + add_dim
+    return Nx_comp, Ny_comp, Nt_comp
+
+
+def generate_coarse_noise(prng, params, data_type, add_dim):
+    """
+    This function generates a spatially and temporally correlated noise.
+    Because it may take a lot of time to compute a correlated noise on
+    a too fine mesh, we recommend to first compute a correlated signal
+    on a coarse mesh, and then use the interpolation function.
+
+    Input:
+        params: (dict) Defines the mesh dimensions and
+            precision. Also define the correlation scales
+
+    Output:
+        (np.array) 3D autocorrelated noise
+    """
 
+    Nx_comp, Ny_comp, Nt_comp = get_nx_ny_nt(data_type, params, add_dim)
+    
     # Generate gaussian noise input·
     #output = np.random.normal(0, 1, (Nx_comp, Ny_comp, Nt_comp))
     output = prng.normal(0, 1, (Nx_comp, Ny_comp, Nt_comp))
diff --git a/chronix2grid/generation/renewable/RenewableBackend.py b/chronix2grid/generation/renewable/RenewableBackend.py
index e7f1a93..35cdd2f 100644
--- a/chronix2grid/generation/renewable/RenewableBackend.py
+++ b/chronix2grid/generation/renewable/RenewableBackend.py
@@ -38,10 +38,21 @@ def __init__(self, out_path, seed, params, loads_charac, res_config_manager, wri
         self.seed = seed
         self.out_path = out_path
 
-    def run(self, solar_pattern=None):
+    def run(self,
+            solar_pattern=None,
+            return_ref_curve=False,
+            return_prng=False,
+            tol_zero=0.):
         """
         Runs the generation model in ``chronix2grid.generation.renewable.generate_solar_wind`` and writes chronics
+        
+        
+        # NB `tol_zero` is set to 0.1 for legacy behaviour, otherwise tests do not pass, but this is a TERRIBLE idea.
         """
         if solar_pattern is None:
             solar_pattern = self.res_config_manager.read_specific()
-        return main(self.out_path, self.seed, self.params, self.loads_charac, solar_pattern, self.write_results)
+        return main(self.out_path, self.seed, self.params, self.loads_charac,
+                    solar_pattern, self.write_results,
+                    return_ref_curve=return_ref_curve,
+                    return_prng=return_prng,
+                    tol_zero=tol_zero)
diff --git a/chronix2grid/generation/renewable/generate_solar_wind.py b/chronix2grid/generation/renewable/generate_solar_wind.py
index 2f9df71..d4a18d3 100755
--- a/chronix2grid/generation/renewable/generate_solar_wind.py
+++ b/chronix2grid/generation/renewable/generate_solar_wind.py
@@ -20,10 +20,30 @@
 import chronix2grid.constants as cst
 
 
-def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, write_results = True):
+def get_add_dim(params, prods_charac):
+    scale_solar_coord_for_correlation = float(params["scale_solar_coord_for_correlation"]) if "scale_solar_coord_for_correlation" in params else None
+    add_dim = 0
+    dx_corr = int(params['dx_corr'])
+    dy_corr = int(params['dy_corr'])
+    for x, y, type_gen  in zip(prods_charac["x"], prods_charac["y"], prods_charac["type"]):
+        if type_gen == "solar" and scale_solar_coord_for_correlation is not None:
+            x = scale_solar_coord_for_correlation * x
+            y = scale_solar_coord_for_correlation * y
+        x_plus = int(x // dx_corr + 1)
+        y_plus = int(y // dy_corr + 1)
+        add_dim = max(y_plus, add_dim)
+        add_dim = max(x_plus, add_dim)
+    return add_dim
+
+
+def main(scenario_destination_path, seed, params, prods_charac, solar_pattern,
+         write_results=True, return_ref_curve=False, return_prng=False,
+         tol_zero=0.):
     """
     This is the solar and wind production generation function, it allows you to generate consumption chronics based on
     production nodes characteristics and on a solar typical yearly production patterns.
+    
+    # NB tol is set to 0.1 for legacy behaviour, otherwise tests do not pass, but this is a TERRIBLE idea.
 
     Parameters
     ----------
@@ -34,13 +54,17 @@ def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, w
     solar_pattern (pandas.DataFrame): hourly solar production pattern for a year. It represent specificity of the production region considered
     smoothdist (float): parameter for smoothing
     write_results (boolean): whether to write results or not. Default is True
-
+    return_ref_curve (bool): whether to return the reference curve or not.
+    return_prng (bool): whether to return (or not) the pseudo random number generator used.
+    
     Returns
     -------
     pandas.DataFrame: solar production chronics generated at every node with additional gaussian noise
     pandas.DataFrame: solar production chronics forecasted for the scenario without additional gaussian noise
     pandas.DataFrame: wind production chronics generated at every node with additional gaussian noise
     pandas.DataFrame: wind production chronics forecasted for the scenario without additional gaussian noise
+
+    And a few other arguments depending on the flags (*eg* return_ref_curve or return_prng)
     """
 
     prng = default_rng(seed)
@@ -69,17 +93,7 @@ def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, w
     # Generate GLOBAL temperature noise
     print('Computing global auto-correlated spatio-temporal noise for sun and wind...')
     scale_solar_coord_for_correlation = float(params["scale_solar_coord_for_correlation"]) if "scale_solar_coord_for_correlation" in params else None
-    add_dim = 0
-    dx_corr = int(params['dx_corr'])
-    dy_corr = int(params['dy_corr'])
-    for x, y, type_gen  in zip(prods_charac["x"], prods_charac["y"], prods_charac["type"]):
-        if type_gen == "solar" and scale_solar_coord_for_correlation is not None:
-            x = scale_solar_coord_for_correlation * x
-            y = scale_solar_coord_for_correlation * y
-        x_plus = int(x // dx_corr + 1)
-        y_plus = int(y // dy_corr + 1)
-        add_dim = max(y_plus, add_dim)
-        add_dim = max(x_plus, add_dim)
+    add_dim = get_add_dim(params, prods_charac)
 
     solar_noise = utils.generate_coarse_noise(prng, params, 'solar', add_dim=add_dim)
     long_scale_wind_noise = utils.generate_coarse_noise(prng, params, 'long_wind', add_dim=add_dim)
@@ -89,12 +103,16 @@ def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, w
     # Compute Wind and solar series of scenario
     print('Generating solar and wind production chronics')
     prods_series = {}
+    solar_ref = None
+    wind_ref = None
+    i_solar = 0
+    i_wind = 0
     for name in prods_charac['name']:
         mask = (prods_charac['name'] == name)
         if prods_charac[mask]['type'].values == 'solar':
             locations = [prods_charac[mask]['x'].values[0], prods_charac[mask]['y'].values[0]]
             Pmax = prods_charac[mask]['Pmax'].values[0]
-            prods_series[name] = swutils.compute_solar_series(
+            tmp_ = swutils.compute_solar_series(
                 prng,
                 locations,
                 Pmax,
@@ -102,12 +120,20 @@ def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, w
                 params, solar_pattern, smoothdist,
                 time_scale=params['solar_corr'],
                 add_dim=add_dim,
-                scale_solar_coord_for_correlation=scale_solar_coord_for_correlation)
-
+                scale_solar_coord_for_correlation=scale_solar_coord_for_correlation,
+                return_ref_curve=return_ref_curve,
+                tol=tol_zero)
+            if return_ref_curve:
+                if solar_ref is None:
+                    solar_ref = np.zeros((tmp_[0].shape[0], (prods_charac['type'].values == 'solar').sum()))
+                prods_series[name], solar_ref[:,i_solar] = tmp_
+            else:
+                prods_series[name] = tmp_
+            i_solar += 1
         elif prods_charac[mask]['type'].values == 'wind':
             locations = [prods_charac[mask]['x'].values[0], prods_charac[mask]['y'].values[0]]
             Pmax = prods_charac[mask]['Pmax'].values[0]
-            prods_series[name] = swutils.compute_wind_series(
+            tmp_ = swutils.compute_wind_series(
                 prng,
                 locations,
                 Pmax,
@@ -115,7 +141,16 @@ def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, w
                 medium_scale_wind_noise,
                 short_scale_wind_noise,
                 params, smoothdist,
-                add_dim=add_dim)
+                add_dim=add_dim,
+                return_ref_curve=return_ref_curve,
+                tol=tol_zero)
+            if return_ref_curve:
+                if wind_ref is None:
+                    wind_ref = np.zeros((tmp_[0].shape[0], (prods_charac['type'].values == 'wind').sum()))
+                prods_series[name], wind_ref[:,i_wind] = tmp_
+            else:
+                prods_series[name] = tmp_
+            i_wind += 1
 
     # Séparation ds séries solaires et éoliennes
     solar_series = {}
@@ -196,5 +231,13 @@ def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, w
             index=False,
             float_format=cst.FLOATING_POINT_PRECISION_FORMAT
         )
-
-    return prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted
\ No newline at end of file
+    if not return_ref_curve and not return_prng:
+        res = prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted
+    elif not return_ref_curve and return_prng:
+        res = prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted, prng
+    elif return_ref_curve and not return_prng:
+        res = prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted, solar_ref, wind_ref
+    elif return_ref_curve and return_prng:
+        res = prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted, solar_ref, wind_ref, prng
+    return res
+    
\ No newline at end of file
diff --git a/chronix2grid/generation/renewable/solar_wind_utils.py b/chronix2grid/generation/renewable/solar_wind_utils.py
index e9c9d0d..60e8b61 100755
--- a/chronix2grid/generation/renewable/solar_wind_utils.py
+++ b/chronix2grid/generation/renewable/solar_wind_utils.py
@@ -15,7 +15,12 @@
 from .. import generation_utils as utils
 import chronix2grid.constants as cst
 
-def compute_wind_series(prng, locations, Pmax, long_noise, medium_noise, short_noise, params, smoothdist, add_dim):
+def compute_wind_series(prng, locations, Pmax, long_noise, medium_noise,
+                        short_noise, params, smoothdist, add_dim,
+                        return_ref_curve=False,
+                        tol=0.):
+    # NB tol is set to 0.0 for legacy behaviour, otherwise tests do not pass, but this is a TERRIBLE idea.
+    
     # Compute refined signals
     long_scale_signal = utils.interpolate_noise(
         long_noise,
@@ -40,7 +45,7 @@ def compute_wind_series(prng, locations, Pmax, long_noise, medium_noise, short_n
     Nt_inter = int(params['T'] // params['dt'] + 1)
     t = np.linspace(0, params['T'], Nt_inter, endpoint=True)
     start_min = int(
-        pd.Timedelta(params['start_date'] - pd.to_datetime('2018/01/01', format='%Y-%m-%d')).total_seconds() // 60)
+        pd.Timedelta(params['start_date'] - pd.to_datetime('2018-01-01', format='%Y-%m-%d')).total_seconds() // 60)
     seasonal_pattern = np.cos((2 * np.pi / (365 * 24 * 60)) * (t - 30 * 24 * 60 - start_min))
 
     # Combine signals
@@ -55,13 +60,21 @@ def compute_wind_series(prng, locations, Pmax, long_noise, medium_noise, short_n
     #signal += np.random.uniform(0, smoothdist, signal.shape) #older version - to be removed
 
     # signal *= 0.95
-    signal[signal < 0.] = 0.
+    signal[signal < tol] = 0.
     signal = smooth(signal)
     wind_series = Pmax * signal
     wind_series[wind_series > 0.95 * Pmax] = 0.95 * Pmax
-    return wind_series
+    if not return_ref_curve:
+        return wind_series
+    return wind_series, 1e-1 * np.exp(4 * (0.7 + 0.3 * seasonal_pattern) * 0.3)
 
-def compute_solar_series(prng, locations, Pmax, solar_noise, params, solar_pattern, smoothdist, time_scale, add_dim, scale_solar_coord_for_correlation=None):
+
+def compute_solar_series(prng, locations, Pmax, solar_noise, params,
+                         solar_pattern, smoothdist, time_scale,
+                         add_dim, scale_solar_coord_for_correlation=None,
+                         return_ref_curve=False,
+                         tol=0.):
+    # NB tol is set to 0.0 for legacy behaviour, otherwise tests do not pass, but this is a TERRIBLE idea.
 
     # Compute noise at desired locations
     if scale_solar_coord_for_correlation is not None:
@@ -69,7 +82,7 @@ def compute_solar_series(prng, locations, Pmax, solar_noise, params, solar_patte
     final_noise = utils.interpolate_noise(solar_noise, params, locations, time_scale, add_dim=add_dim)
 
     # Compute solar pattern
-    solar_pattern = compute_solar_pattern(params, solar_pattern)
+    solar_pattern = compute_solar_pattern(params, solar_pattern, tol=tol)
 
     # Compute solar time series
     std_solar_noise = float(params['std_solar_noise'])
@@ -83,18 +96,23 @@ def compute_solar_series(prng, locations, Pmax, solar_noise, params, solar_patte
     #signal += prng.uniform(0, smoothdist/Pmax, signal.shape) #to be revised: since smmothdist/PMax is very small, the added noise compared to the previous sinal was unsignificant
     #signal += np.random.uniform(0, smoothdist / Pmax, signal.shape) #older version - to be removed
     # signal[signal > 1] = 1
-    signal[signal < 0.] = 0.
+    signal[signal < tol] = 0.
     signal = smooth(signal)
     solar_series = Pmax * signal
     # solar_series[np.isclose(solar_series, 0.)] = 0
     solar_series[solar_series > 0.95 * Pmax] = 0.95 * Pmax
-    return solar_series
+    if not return_ref_curve:
+        return solar_series
+    else:
+        return solar_series, solar_pattern * mean_solar_pattern
 
-def compute_solar_pattern(params, solar_pattern):
+def compute_solar_pattern(params, solar_pattern, tol=0.0):
     """
     Loads a typical hourly pattern, and interpolates it to generate
     a smooth solar generation pattern between 0 and 1
 
+    # NB tol is set to 0.0 for legacy behaviour, otherwise tests do not pass, but this is a TERRIBLE idea.
+    
     Input:
         computation_params: (dict) Defines the mesh dimensions and
             precision. Also define the correlation scales
@@ -104,7 +122,7 @@ def compute_solar_pattern(params, solar_pattern):
         (np.array) A smooth solar pattern
     """
 
-    start_year = pd.to_datetime(str(params['start_date'].year) + '/01/01', format='%Y-%m-%d')
+    start_year = pd.to_datetime(str(params['start_date'].year) + '-01-01', format='%Y-%m-%d')
     end_min = int(pd.Timedelta(params['end_date'] - start_year).total_seconds() // 60)
 
     Nt_inter_hr = int(end_min // 60 + 1)
@@ -112,20 +130,71 @@ def compute_solar_pattern(params, solar_pattern):
     stacked_solar_pattern = solar_pattern
     for i in range(N_repet - 1):
         stacked_solar_pattern = np.append(stacked_solar_pattern, solar_pattern)
-
+    # solar_pattern : over the all year, hourly based
+    
     # The time is in minutes
     t_pattern = 60 * np.linspace(0, 8760 * N_repet, 8760 * N_repet, endpoint=False)
     f2 = interp1d(t_pattern, stacked_solar_pattern, kind='cubic')
+    # f2 = interp1d(t_pattern, stacked_solar_pattern, kind='linear')
 
     Nt_inter = int(params['T'] // params['dt'] + 1)
-    start_year = pd.to_datetime(str(params['start_date'].year) + '/01/01', format='%Y-%m-%d')
+    start_year = pd.to_datetime(str(params['start_date'].year) + '-01-01', format='%Y-%m-%d')
     start_min = int(pd.Timedelta(params['start_date'] - start_year).total_seconds() // 60)
     end_min = int(pd.Timedelta(params['end_date'] - start_year).total_seconds() // 60)
 
     t_inter = np.linspace(start_min, end_min, Nt_inter, endpoint=True)
     output = f2(t_inter)
     output = output * (output > 0)
-
+    output[output < tol] = 0.
+    if "solar_night_hour" in params:
+        # addition: force the solar at 0 at night
+        dts = [params['start_date'] + i * pd.Timedelta(minutes=params['dt']) for i in range(t_inter.shape[0])]
+        dts_hours = np.array([el.hour for el in dts])
+        # i will multiply by 1 almost all hours
+        to_mult = np.ones(len(dts_hours))
+        # except hours provided in solar_night_hour
+        to_mult[np.isin(dts_hours, params["solar_night_hour"])] = 0.
+        
+        # now I want to "smooth" the "to_mult" for the "border"
+        
+        # left smoothing, for 1 day
+        left_ = 1. / (1. + np.exp(1 - ((1+ np.arange(60 // params['dt']) )- (30. /  params['dt']))))
+        # right smoothing (for 1 day too)
+        right_ = 1. / (1. + np.exp(1 + ((1+ np.arange(60 // params['dt']) )- (30. /  params['dt']))))
+        
+        # define some constants
+        nb_day = int((params['end_date'] - params['start_date']).total_seconds() // 86400)
+        solar_night_hour = np.array(params["solar_night_hour"])
+        
+        # handles the "left" border (eg 5 am)
+        first_non_zero_h = np.max(solar_night_hour[solar_night_hour <= 12]) + 1 # 12 for midday
+        all_left = np.tile(left_, nb_day)
+        tmp_ = np.zeros((dts_hours == first_non_zero_h).sum())
+        tmp_[:all_left.shape[0]] = all_left
+        to_mult[dts_hours == first_non_zero_h] = tmp_
+        
+        # handles "right" border (eg 22 / 10 pm)
+        last_non_zero_h = np.min(solar_night_hour[solar_night_hour >= 12]) - 1
+        all_right = np.tile(right_, nb_day)
+        tmp_ = np.zeros((dts_hours == last_non_zero_h).sum())
+        to_rem = tmp_.shape[0] - all_right.shape[0]
+        tmp_[to_rem:] = all_right
+        to_mult[dts_hours == last_non_zero_h] = tmp_
+        
+        # force the output to be 0 at night
+        output *= to_mult
+        
+    if "force_solar_zero" in params:
+        # another heuristic to cap the solar at 0.0 at "night"
+        # when the solar is at 5% of its max value, then 1h after it it should be 0.
+        threshold_zero = output.max() * 0.05
+        max_hour_seen_non_zero = np.max(dts_hours[output >= threshold_zero])
+        max_hour_seen_non_zero += int(params["force_solar_zero"])
+        output[dts_hours >= max_hour_seen_non_zero] = 0.
+        # same in the other direction: 1h before it's always bellow threshold, it should be 0.
+        min_hour_seen_non_zero = np.min(dts_hours[output >= threshold_zero])
+        min_hour_seen_non_zero -= int(params["force_solar_zero"])
+        output[dts_hours <= min_hour_seen_non_zero] = 0.
     return output
 
 def smooth(x, alpha=0.5, beta=None):
diff --git a/chronix2grid/grid2op_utils/add_data.py b/chronix2grid/grid2op_utils/add_data.py
index 849718f..1a626df 100644
--- a/chronix2grid/grid2op_utils/add_data.py
+++ b/chronix2grid/grid2op_utils/add_data.py
@@ -6,7 +6,6 @@
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of Chronix2Grid, A python package to generate "en-masse" chronics for loads and productions (thermal, renewable)
 import os
-import numpy as np
 import json
 from multiprocessing import Pool
 
@@ -14,13 +13,42 @@
 from chronix2grid.grid2op_utils.utils import generate_a_scenario, get_last_scenario_id
 from numpy.random import default_rng
 
-import pdb
+
+# wrapper function for generate_a_scenario
+def generate_a_scenario_wrapper(args):
+    (path_env, name_gen, gen_type, output_dir,
+        start_date, dt, scen_id, load_seed, renew_seed,
+        gen_p_forecast_seed, handle_loss, files_to_copy,
+        save_ref_curve, day_lag, tol_zero, debug) = args
+    res_gen = generate_a_scenario(path_env,
+                                name_gen, gen_type,
+                                output_dir,
+                                start_date, dt,
+                                scen_id,
+                                load_seed, renew_seed, 
+                                gen_p_forecast_seed,
+                                handle_loss,
+                                files_to_copy=files_to_copy,
+                                save_ref_curve=save_ref_curve,
+                                day_lag=day_lag,
+                                tol_zero=tol_zero,
+                                debug=debug)
+    return res_gen
+
 
 def add_data(env: grid2op.Environment.Environment,
              seed=None,
              nb_scenario=1,
              nb_core=1,  # TODO
-             with_loss=True):
+             with_loss=True,
+             files_to_copy=("maintenance_meta.json",
+                            "params_load.json",
+                            "params_forecasts.json"),
+             save_ref_curve=False,
+             day_lag=6,  # TODO 6 because it's 2050
+             debug=False,
+             tol_zero=1e-3
+             ):
     """This function adds some data to already existing scenarios.
     
     .. warning::
@@ -73,7 +101,6 @@ def add_data(env: grid2op.Environment.Environment,
             gen_p_forecast_seeds.append(gen_p_forecast_seed)
     
     # generate the data
-    # TODO multi proc
     path_env = env.get_path_env()
     name_gen = env.name_gen
     gen_type = env.gen_type
@@ -81,6 +108,8 @@ def add_data(env: grid2op.Environment.Environment,
     argss = []
     for j, scen_id in enumerate(scen_ids):
         for i, start_date in enumerate(li_months):
+            # if start_date != "2035-03-12":
+            #     continue
             seed_num = i + j * len(li_months)
             argss.append((path_env,
                           name_gen,
@@ -92,14 +121,16 @@ def add_data(env: grid2op.Environment.Environment,
                           load_seeds[seed_num],
                           renew_seeds[seed_num],
                           gen_p_forecast_seeds[seed_num],
-                          with_loss
+                          with_loss,
+                          files_to_copy,
+                          save_ref_curve,
+                          day_lag,
+                          tol_zero,
+                          debug
                           ))
     if nb_core == 1:
         for args in argss:
-            path_env, name_gen, gen_type, output_dir, start_date, dt, scen_id, load_seed, renew_seed, \
-                gen_p_forecast_seed, handle_loss = args
-            res_gen = generate_a_scenario(path_env, name_gen, gen_type, output_dir, start_date, dt, scen_id, load_seed, renew_seed, 
-                                          gen_p_forecast_seed, handle_loss)
+            res_gen = generate_a_scenario_wrapper(args)
             error_, *_ = res_gen
             if error_ is not None:
                 print("=============================")
@@ -121,4 +152,4 @@ def add_data(env: grid2op.Environment.Environment,
                     json.dump(errors, fp=f)
     else:
         with Pool(nb_core) as p:
-            p.starmap(generate_a_scenario, argss)
+            p.map(generate_a_scenario_wrapper, argss)
diff --git a/chronix2grid/grid2op_utils/gen_utils.py b/chronix2grid/grid2op_utils/gen_utils.py
new file mode 100644
index 0000000..126cbc0
--- /dev/null
+++ b/chronix2grid/grid2op_utils/gen_utils.py
@@ -0,0 +1,408 @@
+# Copyright (c) 2019-2022, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of Chronix2Grid, A python package to generate "en-masse" chronics for loads and productions (thermal, renewable)
+
+import json
+import pandas as pd
+import os
+import cvxpy as cp
+import numpy as np
+
+from chronix2grid.generation.renewable.generate_solar_wind import get_add_dim as get_add_dim_renew
+from chronix2grid.grid2op_utils.noise_generation_utils import (get_forecast,
+                                                               generate_noise)
+
+
+def fix_nan_hydro_i_dont_know_why(arr):
+    # economic_dispatch.make_hydro_constraints_from_res_load_scenario() sometimes return all nans for a given row...
+    # don't ask me why
+    if np.all(np.isfinite(arr)):
+        # nothing to do, data are correct
+        return 1.0 * arr
+    if np.all(~np.isfinite(arr)):
+        # cannot do anything: nothing is finite
+        return arr
+    if np.any(np.isfinite(arr).sum(axis=0) == arr.shape[0]):
+        # there is at least a column full of Nan I cannot do anything either
+        return np.NaN * arr
+    # to fix the nan values, I assign the last known finite value per column
+    res = 1.0 * arr
+    for col_id in range(arr.shape[1]):
+        finite_ind = np.where(np.isfinite(res[:,col_id]))[0]
+        infinite_ind = np.where(~np.isfinite(res[:,col_id]))[0]
+        idx = np.searchsorted(finite_ind, infinite_ind)
+        idx[idx == finite_ind.shape[0]] = finite_ind.shape[0] - 1
+        res[infinite_ind, col_id] = 1.0 * arr[finite_ind[idx], col_id]
+    return res
+
+
+def get_gen_ids_hydro(env):
+    gen_id = 0
+    ids_hyrdo = []
+    for i in range(env.n_gen):
+        if env.gen_redispatchable[i]:
+            if env.gen_type[i] == "hydro":
+                ids_hyrdo.append(gen_id)
+            gen_id += 1
+    ids_hyrdo = np.array(ids_hyrdo)
+    return ids_hyrdo
+    
+
+def fix_forecast_ramps(nb_h,
+                       load_p,
+                       load_p_forecasted,
+                       res_gen_p_df,
+                       res_gen_p_forecasted_df,
+                       env_for_loss,
+                       hydro_constraints,
+                       params):
+    
+    #### cvxpy
+    total_step = 1 # for now
+    total_gen = np.sum(env_for_loss.gen_redispatchable)
+    load_f = load_p_forecasted.sum(axis=1)
+    scaling_factor = env_for_loss.gen_pmax[env_for_loss.gen_redispatchable]
+    scale_for_loads =  np.repeat(scaling_factor.reshape(1,-1), total_step, axis=0)
+    ids_hyrdo = get_gen_ids_hydro(env_for_loss)
+    
+    p_min = np.repeat(env_for_loss.gen_pmin[env_for_loss.gen_redispatchable].reshape(1,-1)  / scaling_factor,
+                      total_step + nb_h,
+                      axis=0)
+    p_max = np.repeat(env_for_loss.gen_pmax[env_for_loss.gen_redispatchable].reshape(1,-1) * params["PmaxErrorCorrRatio"] / scaling_factor,
+                      total_step + nb_h,
+                      axis=0)
+    
+    ramp_min = np.repeat(-env_for_loss.gen_max_ramp_down[env_for_loss.gen_redispatchable].reshape(1,-1) * params["RampErrorCorrRatio"] / scaling_factor,
+                         total_step + (nb_h-1),
+                         axis=0)
+    ramp_max = np.repeat(env_for_loss.gen_max_ramp_up[env_for_loss.gen_redispatchable].reshape(1,-1) * params["RampErrorCorrRatio"] / scaling_factor,
+                         total_step + (nb_h-1),
+                         axis=0)
+    
+    if "hydro_ramp_reduction_factor" in params:
+        ramp_max[:, ids_hyrdo] /= float(params["hydro_ramp_reduction_factor"])
+        ramp_min[:, ids_hyrdo] /= float(params["hydro_ramp_reduction_factor"])
+        
+    res_gen_p = np.zeros((res_gen_p_forecasted_df.shape[0], total_gen))
+    amount_curtailed_for = np.zeros((res_gen_p_forecasted_df.shape[0], ))
+    has_error = np.zeros(res_gen_p_forecasted_df.shape[0], dtype=bool)
+    t0_errors = []
+    errors = []
+    indx_forecasts_for_hydro_only = np.arange(0, nb_h + 1)
+    for t0 in range(res_gen_p_df.shape[0] - 1):
+        # forecast are "consistent from a power system point of view" batch by batch
+        # losses are not handled here !
+        loss_scale = res_gen_p_df.iloc[t0].sum() / load_p.iloc[t0].sum()
+        indx_forecasts = np.arange(t0*nb_h, (t0+1)*nb_h)
+        
+        prod_renew_for = res_gen_p_forecasted_df.iloc[indx_forecasts, env_for_loss.gen_renewable].sum(axis=1)
+        loss_for = (load_f.iloc[indx_forecasts] * loss_scale - prod_renew_for)
+        net_load = np.array(([res_gen_p_df.iloc[t0, env_for_loss.gen_redispatchable].sum()] +  # value in the env
+                             loss_for.tolist())
+                           )
+        # curtailment
+        curt_t = cp.Variable(shape=(total_step + nb_h, ), nonneg=True)
+        # curt_t = cp.multiply(curt_t_scaled, sum_pmax_renew)
+        renew = np.array(([res_gen_p_df.iloc[t0, env_for_loss.gen_renewable].sum()] +  # value in the env
+                          prod_renew_for.tolist())
+                        )
+        scale_curt_factor = np.maximum(renew, 1.)
+        curt_t_scaled = cp.multiply(curt_t, 1. / scale_curt_factor)
+        
+        # generation
+        p_t = cp.Variable(shape=(total_step + nb_h, total_gen), pos=True)
+        real_p = cp.multiply(p_t, scale_for_loads)
+        target_vector = 1.0 * np.concatenate((res_gen_p_df.iloc[t0, env_for_loss.gen_redispatchable].values.reshape(1, total_gen),
+                                              res_gen_p_forecasted_df.iloc[indx_forecasts, env_for_loss.gen_redispatchable].values),
+                                              axis=0).reshape(1+nb_h, total_gen)
+        
+        p_max_here = 1.0 * p_max
+        if hydro_constraints is not None:
+            p_max_here[:, ids_hyrdo] = fix_nan_hydro_i_dont_know_why(hydro_constraints["p_max_pu"].values[indx_forecasts_for_hydro_only, :])
+            indx_forecasts_for_hydro_only += 1
+            indx_forecasts_for_hydro_only[indx_forecasts_for_hydro_only > (res_gen_p_df.shape[0] - 1)] = res_gen_p_df.shape[0] - 1
+        
+        turned_off_orig = 1.0 * (target_vector == 0.)
+        target_vector /= scale_for_loads
+        constraints = [real_p[0,:] == res_gen_p_df.iloc[t0, env_for_loss.gen_redispatchable].values,
+                       p_t >= p_min,
+                       p_t <= p_max_here,
+                       p_t[1:,:] - p_t[:-1,:] >= ramp_min,
+                       p_t[1:,:] - p_t[:-1,:] <= ramp_max,
+                       cp.sum(real_p, axis=1) >= (net_load.reshape(-1) + curt_t),
+                       cp.sum(real_p, axis=1) <= 1.01 * (net_load.reshape(-1) + curt_t),
+                       curt_t <= renew,
+                       curt_t >= 0.,
+                       curt_t[0] == 0.
+                       ]
+        
+        cost = cp.sum_squares(p_t - target_vector) + cp.norm1(cp.multiply(p_t, turned_off_orig)) + 10. * cp.sum_squares(curt_t_scaled)  # TODO normalize last stuff
+        prob = cp.Problem(cp.Minimize(cost), constraints)
+        try:
+            res_opt = prob.solve()
+        except cp.error.SolverError as exc_:
+            t0_errors.append(t0)
+            errors.append(RuntimeError(f"cvxpy failed to find a solution for t0 {t0}, error {exc_}"))
+            has_error[indx_forecasts] = True
+            continue
+        
+        if not np.isfinite(res_opt):
+            t0_errors.append(t0)
+            errors.append(RuntimeError(f"cvxpy failed to find a solution for t0 {t0} and returned an infinite cost"))
+            has_error[indx_forecasts] = True
+            continue
+        
+        # assign the generators
+        gen_p_after_optim = real_p.value
+        if gen_p_after_optim is None:
+            t0_errors.append(t0)
+            errors.append(f"cvxpy failed to find a solution for t0 {t0}, and returned None")
+            has_error[indx_forecasts] = True
+            continue
+        
+        if not has_error[indx_forecasts[0]]:
+            res_gen_p[indx_forecasts] = 1.0 * gen_p_after_optim[1:,:]
+            amount_curtailed_for[indx_forecasts] = curt_t.value[1:]
+            
+    # last value is not used anyway
+    # res_gen_p[-1, :] = 1.0 * gen_p_after_optim[1,:]
+    has_error[-nb_h:] = True
+        
+    res_gen_p_forecasted_df_res = 1.0 * res_gen_p_forecasted_df
+    # assign controlable generators
+    res_gen_p_forecasted_df_res.iloc[~has_error, env_for_loss.gen_redispatchable] = res_gen_p[~has_error,:]
+    
+    # fix for curtailment
+    total_renew = res_gen_p_forecasted_df_res.iloc[:, env_for_loss.gen_renewable].sum(axis=1)
+    total_renew[total_renew <= 1.] = 1.  # normalize for value close to 0.
+    ratio = (total_renew - amount_curtailed_for) / total_renew
+    res_gen_p_forecasted_df_res.iloc[:, env_for_loss.gen_renewable] *= np.stack([ratio for _ in 
+                                                                                 range(np.sum(env_for_loss.gen_renewable))],
+                                                                                 axis=1)
+    
+    # fix for tiny negative value that, with the magic of rounding
+    # can lead to -0.1 in the data
+    res_gen_p_forecasted_df_res[res_gen_p_forecasted_df_res < 0.] = 0.
+    
+    # make sure the forecasts are always above the demands, even the the opf failed
+    return res_gen_p_forecasted_df_res, t0_errors, errors, amount_curtailed_for
+
+
+def fix_negative(forecast):
+    # inplace
+    forecast[forecast < 0.] = 0.
+    
+
+def fix_pmax(forecast, pmax):
+    # forecasts : (nb_gen, nb_t, nb_h)
+    # pmax: (nb_gen)
+    # not inplace
+    pmax_ = np.stack([pmax for _ in range(forecast.shape[1])], axis=1)
+    # pmax_: (nb_gen, nb_t)
+    pmax_ = np.stack([pmax_ for _ in range(forecast.shape[-1])], axis=2)
+    # pmax_: (nb_gen, nb_t, nb_h)
+    
+    init_shape = forecast.shape
+    forecast_ = forecast.ravel()
+    pmax_ = pmax_.ravel()
+    mask_ = forecast_ > pmax_
+    forecast_[mask_] = pmax_[mask_]
+    forecast_ = forecast_.reshape(init_shape)
+    return forecast_
+    
+    
+def fix_solar(forecast, real_value, pmax, tol=1e-2):
+    # inplace !
+    
+    # if value < 0. then it's 0.
+    fix_negative(forecast)
+    
+    # if above pmax then pmax
+    forecast = fix_pmax(forecast, pmax)
+    
+    # no solar at night
+    real_value_ = np.stack([np.roll(real_value, -h) for h in range(forecast.shape[-1])], axis=2)
+    forecast[real_value_ <= tol] = 0.
+
+
+def fix_wind(forecast, real_value, pmax):
+    # inplace !
+    
+    # if value < 0. then it's 0.
+    fix_negative(forecast)
+    
+    # if above pmax then pmax
+    forecast = fix_pmax(forecast, pmax)
+
+
+def generate_new_gen_forecasts(prng,
+                               forecasts_params,
+                               load_params,
+                               loads_charac,
+                               gens_charac,
+                               path_env,
+                               res_gen_p_df):
+    
+    # read the parameters from the inputs
+    nb_gen = len(gens_charac['name'])
+    nb_h = len(forecasts_params["h"])
+    datetime_index = pd.date_range(start=load_params['start_date'],
+                                   end=load_params['end_date'],
+                                   freq=str(load_params['dt']) + 'min')
+    nb_t = datetime_index.shape[0]    
+    
+    # load the parameters controling the RES
+    with open(os.path.join(path_env, "params_res.json"), "r") as f:
+        res_params = json.load(f)
+    
+    res_gen_p_forecasted = None
+    
+    forecasts_params["T"] = load_params["T"]
+    res_params["T"] = load_params["T"]
+    
+    # generate the independant data on the mesh
+    for data_type, fun_fix in zip(["solar", "wind"],
+                                  [fix_solar, fix_wind]):
+        
+        # get the parameters for this data type
+        mask_this_type = gens_charac["type"].values == data_type
+        gen_carac_this_type = gens_charac.loc[mask_this_type]
+        nb_gen_this_type = gen_carac_this_type.shape[0]
+        
+        # generrate the noise for all the generators
+        this_noise, hs, std_hs = generate_noise(loads_charac,
+                                                gens_charac,
+                                                forecasts_params,
+                                                load_params,
+                                                None,
+                                                data_type,
+                                                get_add_dim_renew,
+                                                nb_t,
+                                                res_params,
+                                                gen_carac_this_type,
+                                                nb_gen_this_type,
+                                                add_h0=False,
+                                                prng_noise=prng)
+        # shape: (nb_elem, nb_t, nb_h)
+        
+        # generate all the forecasts
+        gen_p_this_type = 1.0 * res_gen_p_df.loc[:, mask_this_type].values.T
+        gen_p_for_this_type = get_forecast(gen_p_this_type,
+                                           this_noise[:, :-1, :],
+                                           hs,
+                                           std_hs,
+                                           gen_carac_this_type,
+                                           reshape=False,
+                                           keep_first_dim=True)
+        
+        # fix the value of the forecast (if above pmax or bellow pmin for example)
+        fun_fix(gen_p_for_this_type, gen_p_this_type, gen_carac_this_type["Pmax"].values)
+        
+        # now put everything in the right shape
+        if res_gen_p_forecasted is None:
+            res_gen_p_forecasted = np.stack([res_gen_p_df.values.T for _ in hs], axis=2)
+        res_gen_p_forecasted[mask_this_type, :, :] = gen_p_for_this_type
+        
+    res_gen_p_forecasted = res_gen_p_forecasted.reshape(nb_gen, (nb_t - 1) * nb_h)
+    res_gen_p_forecasted = np.transpose(res_gen_p_forecasted, (1, 0))
+    res_gen_p_forecasted_df = pd.DataFrame(res_gen_p_forecasted, columns=gens_charac["name"])
+    return res_gen_p_forecasted_df, nb_h
+
+
+def generate_forecasts_gen(new_forecasts,
+                           prng,
+                           load_p,
+                           load_p_forecasted, 
+                           res_gen_p_df,
+                           sigma,
+                           env_for_loss,
+                           hydro_constraints,
+                           forecasts_params,
+                           load_params,
+                           loads_charac,
+                           gens_charac,
+                           path_env,
+                           opf_params):
+    if new_forecasts:
+        res_gen_p_forecasted_df, nb_h = generate_new_gen_forecasts(prng,
+                                                                   forecasts_params,
+                                                                   load_params,
+                                                                   loads_charac,
+                                                                   gens_charac,
+                                                                   path_env,
+                                                                   res_gen_p_df)
+    else:
+        res_gen_p_forecasted_df = res_gen_p_df * prng.lognormal(mean=0.0,
+                                                                sigma=sigma,
+                                                                size=res_gen_p_df.shape)
+        res_gen_p_forecasted_df = res_gen_p_forecasted_df.shift(-1)
+        res_gen_p_forecasted_df.iloc[-1] = 1.0 * res_gen_p_forecasted_df.iloc[-2]
+        nb_h = 1
+    
+    # "fix" cases where forecasts are bellow the loads => in that case scale the
+    # controlable generation to be at least 1% above total demand
+    total_gen = res_gen_p_forecasted_df.sum(axis=1)
+    total_demand = load_p_forecasted.sum(axis=1)
+    mask_ko = total_gen <= total_demand
+    nb_concerned = (mask_ko).sum()
+    tmp = type(env_for_loss).gen_pmax[env_for_loss.gen_redispatchable]
+    tmp = tmp / tmp.sum()
+    rep_factor = np.tile(tmp.reshape(-1,1), nb_concerned).T
+    res_gen_p_forecasted_df.loc[mask_ko, type(env_for_loss).gen_redispatchable] *= (1.01 * total_demand - total_gen)[mask_ko].values.reshape(-1,1) * rep_factor
+    
+    # and fix the ramps (an optimizer, step by step)
+    tmp_ = fix_forecast_ramps(nb_h,
+                              load_p,
+                              load_p_forecasted,
+                              res_gen_p_df,
+                              res_gen_p_forecasted_df,
+                              env_for_loss,
+                              hydro_constraints,
+                              opf_params)
+    res_gen_p_forecasted_df_res, t0_errors, errors, amount_curtailed_for = tmp_
+    return res_gen_p_forecasted_df_res, amount_curtailed_for, t0_errors, errors
+
+
+def apply_maintenance_wind_farm(extra_winds_params, prod_wind_init,
+                                start_date_dt, end_date_dt, dt_min,
+                                renew_prng):
+    res = 1.0 * prod_wind_init
+    # perfom some checks to be sure
+    assert "proba_daily_outage_per_wind_mill" in extra_winds_params  # proba for the time
+    assert "nb_wind_turbine" in extra_winds_params
+    assert "p_geom_failure" in extra_winds_params # proba for the number
+    assert len(extra_winds_params["nb_wind_turbine"]) == prod_wind_init.shape[1]
+    for el in res.columns:
+        assert el in extra_winds_params["nb_wind_turbine"], f"error {el} is not in extra_winds_params.json `nb_wind_turbine`"
+    nb_windturb_per_mill = [extra_winds_params["nb_wind_turbine"][el] for el in res.columns]
+    proba_wind_turb_fails = float(extra_winds_params["p_geom_failure"])
+    
+    # now simulate the "failure" process
+    proba_failure_per_time = float(extra_winds_params["proba_daily_outage_per_wind_mill"]) / (60. * 24. / float(dt_min))
+    do_fail = renew_prng.uniform(size=res.shape)
+    do_fail = do_fail <= proba_failure_per_time  # whether or not an external event caused a failure
+    nb_failure = do_fail.sum(axis=0)  # number of failure in the scenario
+    mult_factor = np.ones(res.shape)
+    for wind_turb_id in range(res.shape[1]):
+        nb_turb_in_service = np.full(shape=res.shape[0], fill_value=nb_windturb_per_mill[wind_turb_id])
+        time_start_failure = np.where(do_fail[:,wind_turb_id])[0]
+        
+        # simulate intensity and duration of failures
+        fail_id = 0
+        while fail_id < nb_failure[wind_turb_id]:
+            ts_start = time_start_failure[fail_id]
+            # intensity of the external event : number of wind_mill affected
+            do_turbine_fail = renew_prng.uniform(size=nb_turb_in_service[ts_start]) <=  proba_wind_turb_fails
+            # duration of the outage
+            duration = int(60 * 24 / float(dt_min))  # fixed at 1 day for now  # TODO !
+            # and now update everything
+            duration_ids = np.arange(duration) + ts_start
+            duration_ids = duration_ids[duration_ids < res.shape[0]]
+            nb_turb_in_service[duration_ids] -= do_turbine_fail.sum()
+            fail_id += 1
+        mult_factor[:, wind_turb_id] = nb_turb_in_service / nb_windturb_per_mill[wind_turb_id]
+    return res * mult_factor
diff --git a/chronix2grid/grid2op_utils/generate_one_episode.py b/chronix2grid/grid2op_utils/generate_one_episode.py
index 4856a62..0d6e1f7 100644
--- a/chronix2grid/grid2op_utils/generate_one_episode.py
+++ b/chronix2grid/grid2op_utils/generate_one_episode.py
@@ -20,7 +20,8 @@ def generate_one_episode(env: grid2op.Environment.Environment,
                          start_date,
                          nb_steps=None,
                          seed=None,
-                         with_loss=True):
+                         with_loss=True,
+                         files_to_copy=("maintenance_meta.json", "params_load.json", "params_forecasts.json")):
     """This function adds some data to already existing scenarios.
     
     .. warning::
@@ -57,6 +58,7 @@ def generate_one_episode(env: grid2op.Environment.Environment,
     output_dir = None
     while error_ is not None:
         res_gen = generate_a_scenario(path_env, name_gen, gen_type, output_dir, start_date, dt, scen_id, load_seed, renew_seed, 
-                                      gen_p_forecast_seed, with_loss, nb_steps=nb_steps)
+                                      gen_p_forecast_seed, with_loss, nb_steps=nb_steps,
+                                      files_to_copy=files_to_copy)
         error_, quality_, load_p, load_p_forecasted, load_q, load_q_forecasted, res_gen_p_df, res_gen_p_forecasted_df = res_gen
     return load_p, load_p_forecasted, load_q, load_q_forecasted, res_gen_p_df, res_gen_p_forecasted_df
diff --git a/chronix2grid/grid2op_utils/loads_utils.py b/chronix2grid/grid2op_utils/loads_utils.py
new file mode 100644
index 0000000..cff13ae
--- /dev/null
+++ b/chronix2grid/grid2op_utils/loads_utils.py
@@ -0,0 +1,207 @@
+# Copyright (c) 2019-2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of Chronix2Grid, A python package to generate "en-masse" chronics for loads and productions (thermal, renewable)
+
+import os
+import json
+import pandas as pd
+import numpy as np
+
+from chronix2grid.getting_started.example.input.generation.patterns import ref_pattern_path
+from chronix2grid.generation.consumption import ConsumptionGeneratorBackend
+from chronix2grid.generation.consumption.consumption_utils import (get_seasonal_pattern,
+                                                                    compute_load_pattern)
+
+from chronix2grid.generation.consumption.generate_load import get_add_dim as get_add_dim_load
+
+from chronix2grid.grid2op_utils.noise_generation_utils import (get_forecast,
+                                                               generate_noise)
+       
+       
+def get_data_frames(load_p, load_p_for, loads_charac, datetime_index, nb_h):
+    load_p_for = np.transpose(load_p_for, (1, 0))
+    load_p = np.transpose(load_p, (1, 0))
+    load_p_df = pd.DataFrame(load_p, columns=loads_charac["name"])
+    load_p_for_df = pd.DataFrame(load_p_for, columns=loads_charac["name"])
+    
+    load_p_df["datetime"] = datetime_index
+    load_p_df.set_index('datetime', inplace=True)
+    load_p_df = load_p_df.sort_index(ascending=True)
+    size_after = len(load_p_df) - 1
+    load_p_df = load_p_df.head(size_after)
+    
+    load_p_for_df = load_p_for_df.head(size_after * nb_h)
+    return load_p_df, load_p_for_df
+    
+    
+def get_load_ref(loads_charac, load_params, load_weekly_pattern, isoweekday_lwp=None, hour_minutes_lwp=None):
+    weekly_pattern = load_weekly_pattern['test'].values
+    pmax_weekly_pattern = []
+    for index, name in enumerate(loads_charac['name']):
+        mask = (loads_charac['name'] == name)
+        Pmax = loads_charac[mask]['Pmax'].values[0]
+        tmp_ = Pmax * compute_load_pattern(load_params, weekly_pattern, index,
+                                           isoweekday_lwp=isoweekday_lwp,
+                                           hour_minutes_lwp=hour_minutes_lwp)
+        pmax_weekly_pattern.append(tmp_.reshape(1, -1))
+    load_ref = np.concatenate(pmax_weekly_pattern)
+    return load_ref
+
+
+def generate_new_loads(load_seed,
+                       load_params,
+                       forecasts_params,
+                       loads_charac,
+                       gen_charac,
+                       load_weekly_pattern,
+                       data_type='temperature',
+                       day_lag=6,
+                       return_ref_curve=True):    
+    # read the parameters from the inputs
+    nb_load = len(loads_charac['name'])
+    nb_h = len(forecasts_params["h"])
+    datetime_index = pd.date_range(start=load_params['start_date'],
+                                   end=load_params['end_date'],
+                                   freq=str(load_params['dt']) + 'min')
+    nb_t = datetime_index.shape[0] 
+    std_temperature_noise = float(load_params['std_temperature_noise'])
+    
+    # retrieve the reference curve "bar"
+    datetime_lwp = pd.to_datetime(load_weekly_pattern["datetime"], format="%Y-%m-%d %H:%M:%S")
+    isoweekday = np.array([el.isoweekday() for el in datetime_lwp])
+    hour_minutes = np.array([el.hour * 60 + el.minute for el in datetime_lwp])
+    load_ref = get_load_ref(loads_charac, load_params, load_weekly_pattern,
+                            isoweekday_lwp=isoweekday, hour_minutes_lwp=hour_minutes)
+    # (nb_load, nb_t)
+    
+    # Compute seasonal pattern
+    seasonal_pattern_unit = get_seasonal_pattern(load_params)
+    seasonal_pattern_load = np.tile(seasonal_pattern_unit, (nb_load, 1))
+    # (nb_load, nb_t)
+    if return_ref_curve:
+        load_hat = load_ref * seasonal_pattern_load
+        load_hat = 1.0 * load_hat.T[:-1,:]
+    
+    # compute the noise for the loads
+    loads_noise, hs, std_hs = generate_noise(loads_charac,
+                                             gen_charac,
+                                             forecasts_params,
+                                             load_params,
+                                             load_seed,
+                                             data_type,
+                                             get_add_dim_load,
+                                             nb_t,
+                                             load_params,
+                                             loads_charac,
+                                             nb_load,
+                                             add_h0=True)
+        
+    # generate the "real" loads
+    load_p = load_ref * (std_temperature_noise * loads_noise[:,:,0] + seasonal_pattern_load)
+    # shape (n_load, nb_t)
+    
+    # now generate all the forecasts
+    load_p_for = get_forecast(load_p, loads_noise, hs, std_hs, loads_charac)
+    # shape (n_load, nb_t, nb_h)
+    
+    # create the data frames
+    load_p_df, load_p_for_df = get_data_frames(load_p,
+                                               load_p_for,
+                                               loads_charac,
+                                               datetime_index,
+                                               nb_h)
+    if not return_ref_curve:
+        return load_p_df, load_p_for_df
+    return load_p_df, load_p_for_df, load_hat
+
+
+def generate_loads(path_env,
+                   load_seed,
+                   forecast_prng,
+                   start_date_dt,
+                   end_date_dt,
+                   dt,
+                   number_of_minutes,
+                   generic_params,
+                   load_q_from_p_coeff_default=0.7,
+                   day_lag=6):
+    """
+    This function generates the load for each consumption on a grid
+
+    Parameters
+    ----------
+    path_env : _type_
+        _description_
+    load_seed : _type_
+        _description_
+    start_date_dt : _type_
+        _description_
+    end_date_dt : _type_
+        _description_
+    dt : _type_
+        _description_
+    number_of_minutes : _type_
+        _description_
+    generic_params : _type_
+        _description_
+
+    Returns
+    -------
+    _type_
+        _description_
+    """
+    with open(os.path.join(path_env, "params_load.json"), "r") as f:
+        load_params = json.load(f)
+    load_params["start_date"] = start_date_dt
+    load_params["end_date"] = end_date_dt
+    load_params["dt"] = int(dt)
+    load_params["T"] = number_of_minutes
+    load_params["planned_std"] = float(generic_params["planned_std"])
+    
+    if "load_q_from_p_coeff" in load_params:
+        load_q_from_p_coeff = float(load_params["load_q_from_p_coeff"])
+    else:
+        load_q_from_p_coeff = load_q_from_p_coeff_default
+        
+    forecasts_params = {}
+    new_forecasts = False
+    path_for_ = os.path.join(path_env, "params_forecasts.json")
+    if os.path.exists(path_for_):
+        with open(path_for_, "r") as f:
+            forecasts_params = json.load(f)
+        new_forecasts = True
+    
+    loads_charac = pd.read_csv(os.path.join(path_env, "loads_charac.csv"), sep=",")
+    gen_charac = pd.read_csv(os.path.join(path_env, "prods_charac.csv"), sep=",")
+    load_weekly_pattern = pd.read_csv(os.path.join(ref_pattern_path, "load_weekly_pattern.csv"), sep=",")
+    
+    if new_forecasts:
+        load_p, load_p_forecasted, load_ref_curve = generate_new_loads(load_seed,
+                                                                       load_params,
+                                                                       forecasts_params,
+                                                                       loads_charac,
+                                                                       gen_charac,
+                                                                       load_weekly_pattern,
+                                                                       day_lag=day_lag,
+                                                                       return_ref_curve=True)
+    else:
+        load_generator = ConsumptionGeneratorBackend(out_path=None,
+                                                     seed=load_seed, 
+                                                     params=load_params,
+                                                     loads_charac=loads_charac,
+                                                     write_results=False,
+                                                     load_config_manager=None,
+                                                     day_lag=day_lag)
+        
+        load_p, load_p_forecasted, load_ref_curve = load_generator.run(load_weekly_pattern=load_weekly_pattern,
+                                                                       return_ref_curve=True,
+                                                                       use_legacy=False)
+    
+    load_q = load_p * load_q_from_p_coeff
+    load_q_forecasted = load_p_forecasted * load_q_from_p_coeff
+    return (new_forecasts, forecasts_params, load_params, loads_charac,
+            load_p, load_q, load_p_forecasted, load_q_forecasted, load_ref_curve)
diff --git a/chronix2grid/grid2op_utils/noise_generation_utils.py b/chronix2grid/grid2op_utils/noise_generation_utils.py
new file mode 100644
index 0000000..9837a24
--- /dev/null
+++ b/chronix2grid/grid2op_utils/noise_generation_utils.py
@@ -0,0 +1,238 @@
+# Copyright (c) 2019-2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of Chronix2Grid, A python package to generate "en-masse" chronics for loads and productions (thermal, renewable)
+
+import itertools
+import numpy as np
+import warnings
+
+from chronix2grid.generation.generation_utils import get_nx_ny_nt
+from numpy.random import default_rng
+
+
+def generate_coords_mesh(Nx_comp, Ny_comp, Nt_comp, Nh_comp, size,
+                         ratio_adjust=1.0):
+    # distort the distance for the mesh (because we use kNN later on)
+    max_mesh = max(Nx_comp, Ny_comp, Nt_comp, Nh_comp)
+    rho_mesh_x = max_mesh / Nx_comp * ratio_adjust
+    rho_mesh_y = max_mesh / Ny_comp * ratio_adjust
+    rho_mesh_t = max_mesh / Nt_comp * ratio_adjust
+    rho_mesh_h = max_mesh / Nh_comp * ratio_adjust
+    
+    # generate coordinates of the mesh
+    coords_mesh = np.empty((size, 4))
+    for i, row in enumerate(itertools.product(range(Nx_comp), range(Ny_comp), range(Nt_comp), range(Nh_comp))):
+        coords_mesh[i,:] = row
+    coords_mesh[:,0] /= (Nx_comp - 1) * rho_mesh_x
+    coords_mesh[:,1] /= (Ny_comp - 1) * rho_mesh_y
+    coords_mesh[:,2] /= (Nt_comp - 1) * rho_mesh_t
+    coords_mesh[:,3] /= (Nh_comp - 1) * rho_mesh_h
+    return coords_mesh, rho_mesh_x, rho_mesh_y, rho_mesh_t, rho_mesh_h
+    
+
+def get_load_mesh_tmp(nb_t, hs_, hs_mins, rho_mesh_t, rho_mesh_h):
+    nb_h = len(hs_)
+    # now find the noise for all the loads
+    load_mesh_tmp = np.empty((nb_t * nb_h, 4))
+    for i, row in enumerate(itertools.product(range(nb_t), hs_)):
+        load_mesh_tmp[i,:] = (0, 0, *row)
+    # compute where the "point" is on the mesh
+    load_mesh_tmp[:,2] += 1
+    load_mesh_tmp[:,2] /= (nb_t + 1)
+    load_mesh_tmp[:,3] /= (np.max(hs_mins) + 5)
+    
+    # don't forget to re normalize it
+    load_mesh_tmp[:,2] /= rho_mesh_t
+    load_mesh_tmp[:,3] /= rho_mesh_h
+    return load_mesh_tmp
+
+
+def compute_noise(load_mesh_tmp,
+                  loads_charac,
+                  model,
+                  range_x, range_y,
+                  delta_x, delta_y,
+                  rho_mesh_x, rho_mesh_y,
+                  nb_t, nb_h, nb_load):
+    
+    loads_noise = np.zeros((nb_load, nb_t, nb_h))
+    for row_id, (load_id, (load_x, load_y)) in enumerate(loads_charac[["x", 'y']].iterrows()):
+        # compute where the "point" is on the mesh
+        load_mesh = 1.0 * load_mesh_tmp
+        load_mesh[:,0] = (load_x - range_x[0])/delta_x
+        load_mesh[:,1] = (load_y - range_y[0])/delta_y
+        
+        # don't forget to re normalize it
+        load_mesh[:,0] /= rho_mesh_x
+        load_mesh[:,1] /= rho_mesh_y
+        
+        noise_load = model.predict(load_mesh)
+        loads_noise[row_id, :, :] = noise_load.reshape(nb_t, nb_h)
+        
+    # renormalize the noise
+    loads_noise -= loads_noise.mean()
+    loads_noise /= loads_noise.std()  
+    return loads_noise
+
+
+def get_knn_fitted(forecasts_params, coords_mesh, noise_mesh):
+    from sklearn.neighbors import KNeighborsRegressor
+    if "n_neighbors" in forecasts_params:
+        n_neighbors = int(forecasts_params["n_neighbors"])
+        assert n_neighbors > 0, f"n_neighbors parameters should be > 0 and is currently {n_neighbors}"
+    else:
+        n_neighbors = 16 * 15 - 32  # 16 = les 16 arretes du cube 4d dans lequel il est
+    
+    if "leaf_size_knn" in forecasts_params:
+        leaf_size_knn = int(forecasts_params["leaf_size_knn"])
+    else:
+        leaf_size_knn = 100  # found to be the best compromise
+    
+    if "algorithm_knn" in forecasts_params:
+        algorithm_knn = str(forecasts_params["algorithm_knn"])
+    else:
+        algorithm_knn = "auto"  # should be good
+        
+    model = KNeighborsRegressor(n_neighbors=n_neighbors,  
+                                weights="distance",
+                                leaf_size=leaf_size_knn,
+                                algorithm=algorithm_knn)
+    model.fit(coords_mesh, noise_mesh.reshape(-1))
+    return model
+
+
+def get_forecast_parameters(forecasts_params, load_params, data_type="load"):    
+    hs_mins = [int(el) for el in forecasts_params["h"]]
+    # convert the h in "number of steps" (and not in minutes)
+    hs = [h // load_params['dt'] for h in hs_mins]
+    # check consistency
+    for h, h_min in zip(hs, hs_mins):
+        if h * load_params['dt'] != h_min:
+            # TODO logger
+            warnings.warn("Some forecast horizons are not muliple of the duration of a steps. They will be rounded")
+    
+    # load the parameters "h" for the forecasts
+    std_hs = [float(el) for el in forecasts_params[f"h_std_{data_type}"]]
+    for std in std_hs:
+        assert std > 0, f"all parameters of 'h_std_{data_type}' should be >0 we found {std}"
+    assert len(std_hs) == len(hs), "you should provide as much 'error' as there are forecasts horizon"
+    return hs_mins, hs, std_hs
+
+
+def get_iid_noise(load_seed, load_params, forecasts_params,
+                  loads_charac, data_type, get_add_dim,
+                  prng=None):
+    if prng is None:
+        prng = default_rng(load_seed)
+    add_dim = get_add_dim(load_params, loads_charac)
+    Nx_comp, Ny_comp, Nt_comp = get_nx_ny_nt(data_type, load_params, add_dim)
+    Nh_comp = forecasts_params["nb_h_iid"]
+    noise_mesh = prng.normal(0, 1, (Nx_comp, Ny_comp, Nt_comp, Nh_comp))
+    return noise_mesh, (Nx_comp, Ny_comp, Nt_comp, Nh_comp)
+
+
+def resize_mesh_factor(loads_charac, gen_charac, ratio_border=20.):
+    # in the generation we will map this mesh to [0., 1.]
+    range_x = min(loads_charac["x"].min(), gen_charac["x"].min()), max(loads_charac["x"].max(), gen_charac["x"].max())
+    range_y = min(loads_charac["y"].min(), gen_charac["y"].min()), max(loads_charac["y"].max(), gen_charac["y"].max())
+    delta_x_ = range_x[1] - range_x[0]
+    range_x = (range_x[0] - delta_x_ / ratio_border, range_x[1] + delta_x_ / ratio_border)
+    delta_y_ = range_y[1] - range_y[0]
+    range_y = (range_y[0] - delta_y_/ ratio_border, range_y[1] + delta_y_ / ratio_border)
+    # the (+/- delta_x_ / 20.) is to handle side effect at the border of the mesh
+    
+    delta_x = range_x[1] - range_x[0]
+    delta_y = range_y[1] - range_y[0]
+    return delta_x, delta_y, range_x, range_y
+
+    
+def get_forecast(load_p, loads_noise, hs, std_hs, loads_charac, reshape=True, keep_first_dim=False):
+    nb_load = load_p.shape[0]
+    nb_t = load_p.shape[1]
+    nb_h = len(hs)
+    
+    load_p_for = np.stack([np.roll(load_p, -h) for h in hs], axis=2)
+    
+    if keep_first_dim:
+        loads_noise_forecast = 1.0 * loads_noise[:,:,:]
+    else:
+        # in case of load, first dimension is thrown out (used to
+        # compute the real time series data)
+        loads_noise_forecast = 1.0 * loads_noise[:,:,1:]
+    loads_noise_forecast *= np.array(std_hs).reshape(1, len(hs))
+    loads_noise_forecast *= loads_charac["Pmax"].values.reshape(nb_load, 1, 1)
+    load_p_for += loads_noise_forecast
+    # keep the last value for the "forecasts outside the episode"
+    for col, ts in enumerate(hs):
+        load_p_for[:, (nb_t - ts):, col] = load_p_for[:, (nb_t - ts) - 1, col].reshape(nb_load, 1)
+    
+    # put everything in the right shape
+    if reshape:
+        load_p_for = load_p_for.reshape(nb_load, nb_t * nb_h)
+    return load_p_for
+
+
+def generate_noise(loads_charac,
+                   gen_charac,
+                   forecasts_params,
+                   load_params,
+                   load_seed,
+                   data_type,
+                   get_add_dim,
+                   nb_t,
+                   elem_params,
+                   elem_charac,
+                   nb_elem,
+                   add_h0=True,
+                   prng_noise=None):
+    # compute the "real" size of the mesh 
+    delta_x, delta_y, range_x, range_y = resize_mesh_factor(loads_charac, gen_charac)
+    
+    # retrieve the forecast parameters
+    if data_type=='temperature':
+        hs_mins, hs, std_hs = get_forecast_parameters(forecasts_params, load_params) 
+    else:
+        hs_mins, hs, std_hs = get_forecast_parameters(forecasts_params, load_params, data_type)
+             
+    # generate the independant data on the mesh 
+    tmp_ = get_iid_noise(load_seed, elem_params, forecasts_params,
+                         loads_charac, data_type, get_add_dim,
+                         prng_noise)
+    noise_mesh, (Nx_comp, Ny_comp, Nt_comp, Nh_comp) = tmp_
+    
+    # resize the mesh
+    res_mesh = generate_coords_mesh(Nx_comp,
+                                    Ny_comp,
+                                    Nt_comp,
+                                    Nh_comp,
+                                    noise_mesh.size) 
+    coords_mesh, rho_mesh_x, rho_mesh_y, rho_mesh_t, rho_mesh_h = res_mesh
+    
+    # "fit" the kNN    
+    model = get_knn_fitted(forecasts_params, coords_mesh, noise_mesh)
+    
+    # get the "temporary" for the load coordinates (whether or not to generate the
+    # also the "env" value with this function)
+    if add_h0:
+        hs_ = [0] + hs
+    else:
+        hs_ = hs
+    
+    # retrieve the "template" of coordinates of the mesh for the output vector
+    mesh_tmp = get_load_mesh_tmp(nb_t, hs_, hs_mins,
+                                 rho_mesh_t, rho_mesh_h)
+    
+    # now retrieve the real noise for each load
+    this_noise = compute_noise(mesh_tmp,
+                               elem_charac,
+                               model,
+                               range_x, range_y,
+                               delta_x, delta_y,
+                               rho_mesh_x, rho_mesh_y,
+                               nb_t, len(hs_), nb_elem)
+    
+    return this_noise, hs, std_hs
\ No newline at end of file
diff --git a/chronix2grid/grid2op_utils/utils.py b/chronix2grid/grid2op_utils/utils.py
index 39a26df..ea42ca6 100644
--- a/chronix2grid/grid2op_utils/utils.py
+++ b/chronix2grid/grid2op_utils/utils.py
@@ -16,20 +16,25 @@
 import pandas as pd
 import os
 import cvxpy as cp
+import numpy as np
+from numpy.random import default_rng
 import grid2op
 from grid2op.Parameters import Parameters
-from grid2op.Chronics import ChangeNothing, FromNPY
+from grid2op.Chronics import FromNPY
 from grid2op.Action import DontAct
 from grid2op.Opponent import NeverAttackBudget, BaseOpponent
 from lightsim2grid import LightSimBackend
-import numpy as np
-from numpy.random import default_rng
 
-from chronix2grid.generation.consumption import ConsumptionGeneratorBackend
+
 from chronix2grid.generation.renewable import RenewableBackend
 from chronix2grid.generation.dispatch.PypsaDispatchBackend import PypsaDispatcher
 from chronix2grid.getting_started.example.input.generation.patterns import ref_pattern_path
 from chronix2grid.generation.dispatch.EconomicDispatch import ChroniXScenario
+from chronix2grid.grid2op_utils.loads_utils import generate_loads
+from chronix2grid.grid2op_utils.gen_utils import (generate_forecasts_gen,
+                                                  fix_nan_hydro_i_dont_know_why,
+                                                  get_gen_ids_hydro,
+                                                  apply_maintenance_wind_farm)
 
 import warnings
 
@@ -77,62 +82,10 @@ def get_last_scenario_id(env_chronics_dir):
     if max_ is None:
         max_ = -1
     return max_
-            
-    
-def generate_loads(path_env, load_seed, start_date_dt, end_date_dt, dt, number_of_minutes, generic_params,
-                   load_q_from_p_coeff=0.7,
-                   day_lag=6):
-    """
-    This function generates the load for each consumption on a grid
-
-    Parameters
-    ----------
-    path_env : _type_
-        _description_
-    load_seed : _type_
-        _description_
-    start_date_dt : _type_
-        _description_
-    end_date_dt : _type_
-        _description_
-    dt : _type_
-        _description_
-    number_of_minutes : _type_
-        _description_
-    generic_params : _type_
-        _description_
-
-    Returns
-    -------
-    _type_
-        _description_
-    """
-    with open(os.path.join(path_env, "params_load.json"), "r") as f:
-        load_params = json.load(f)
-    load_params["start_date"] = start_date_dt
-    load_params["end_date"] = end_date_dt
-    load_params["dt"] = int(dt)
-    load_params["T"] = number_of_minutes
-    load_params["planned_std"] = float(generic_params["planned_std"])
-    
-    loads_charac = pd.read_csv(os.path.join(path_env, "loads_charac.csv"), sep=",")
-    load_weekly_pattern = pd.read_csv(os.path.join(ref_pattern_path, "load_weekly_pattern.csv"), sep=",")
-    
-    load_generator = ConsumptionGeneratorBackend(out_path=None,
-                                                 seed=load_seed, 
-                                                 params=load_params,
-                                                 loads_charac=loads_charac,
-                                                 write_results=False,
-                                                 load_config_manager=None,
-                                                 day_lag=day_lag)
-    
-    load_p, load_p_forecasted = load_generator.run(load_weekly_pattern=load_weekly_pattern)
-    load_q = load_p * load_q_from_p_coeff
-    load_q_forecasted = load_p_forecasted * load_q_from_p_coeff
-    return load_p, load_q, load_p_forecasted, load_q_forecasted
 
-
-def generate_renewable_energy_sources(path_env, renew_seed, start_date_dt, end_date_dt, dt, number_of_minutes, generic_params, gens_charac):
+def generate_renewable_energy_sources(path_env, renew_seed, start_date_dt, end_date_dt,
+                                      dt, number_of_minutes, generic_params, gens_charac,
+                                      tol_zero=1e-3):
     """This function generates the amount of power produced by renewable energy sources (res). 
     
     It serves as a maximum value for the economic dispatch. 
@@ -175,12 +128,18 @@ def generate_renewable_energy_sources(path_env, renew_seed, start_date_dt, end_d
                                      loads_charac=gens_charac,
                                      res_config_manager=None,
                                      write_results=False)
-    prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted = renew_backend.run(solar_pattern=solar_pattern)
-    return prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted
+    tmp_ = renew_backend.run(solar_pattern=solar_pattern,
+                             return_ref_curve=True,
+                             return_prng=True,
+                             tol_zero=tol_zero)
+    prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted, solar_ref, wind_ref, renew_prng = tmp_
+    return prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted, solar_ref, wind_ref, renew_prng
 
 
 def generate_economic_dispatch(path_env, start_date_dt, end_date_dt, dt, number_of_minutes, generic_params, 
-                               load_p, prod_solar, prod_wind, name_gen, gen_type, scenario_id, final_gen_p, gens_charac):
+                               load_p, prod_solar, prod_wind, name_gen, gen_type, scenario_id, final_gen_p,
+                               gens_charac,
+                               opf_params):
     """This function emulates a perfect market where all productions need to meet the demand at the minimal cost.
     
     It does not consider limit on powerline, nor contigencies etc. The power network does not exist here. Only the ramps and
@@ -220,13 +179,6 @@ def generate_economic_dispatch(path_env, start_date_dt, end_date_dt, dt, number_
     _type_
         _description_
     """
-    with open(os.path.join(path_env, "params_opf.json"), "r") as f:
-        opf_params = json.load(f)
-    opf_params["start_date"] = start_date_dt
-    opf_params["end_date"] = end_date_dt
-    opf_params["dt"] = int(dt)
-    opf_params["T"] = number_of_minutes
-    opf_params["planned_std"] = float(generic_params["planned_std"])
     
     load = pd.DataFrame(load_p.sum(axis=1))
     total_solar = prod_solar.sum(axis=1)
@@ -280,194 +232,9 @@ def generate_economic_dispatch(path_env, start_date_dt, end_date_dt, dt, number_
     
     total_wind_curt = total_wind.values.sum() - res_dispatch.chronix.prods_dispatch['agg_wind'].values.sum()
     total_solar_curt = total_solar.values[mask_solar].sum() - res_dispatch.chronix.prods_dispatch['agg_solar'].values[mask_solar].sum()
-    return final_gen_p, total_wind_curt, total_solar_curt, None
-
-
-def _adjust_gens_old(all_loss_orig,
-                 env_for_loss,
-                 datetimes,
-                 total_solar,
-                 total_wind,
-                 params,
-                 env_path,
-                 env_param,
-                 load_without_loss,
-                 load_p, 
-                 load_q,
-                 gen_p,
-                 gen_v,
-                 economic_dispatch,
-                 diff_,
-                 threshold_stop=0.1,  # stop when all generators move less that this
-                 max_iter=100,  # declare a failure after this number of iteration
-                 iter_quality_decrease=50,  # acept a reduction of the quality after this number of iteration
-                 percentile_quality_decrease=99,
-                 ):
-    """This function is an auxilliary function.
+    return final_gen_p, total_wind_curt, total_solar_curt, hydro_constraints, None
+    
     
-    Like its main one (see handle_losses) it is here to make sure that if you run an AC model with the data generated, 
-    then the generator setpoints will not change too much 
-    (less than `threshold_stop` MW)
-
-    Parameters
-    ----------
-    all_loss_orig : _type_
-        _description_
-    env_for_loss : _type_
-        _description_
-    datetimes : _type_
-        _description_
-    total_solar : _type_
-        _description_
-    total_wind : _type_
-        _description_
-    params : _type_
-        _description_
-    env_path : _type_
-        _description_
-    env_param : _type_
-        _description_
-    load_without_loss : _type_
-        _description_
-    load_p : _type_
-        _description_
-    load_q : _type_
-        _description_
-    gen_p : _type_
-        _description_
-    gen_v : _type_
-        _description_
-    economic_dispatch : _type_
-        _description_
-    diff_ : _type_
-        _description_
-    threshold_stop : float, optional
-        _description_, by default 0.1
-
-    Returns
-    -------
-    _type_
-        _description_
-    """
-    quality_ = None
-    error_ = None
-    if np.any(~np.isfinite(gen_p)):
-        error_ = RuntimeError("Input data contained Nans !")
-        return None, error_, quality_
-    all_loss = all_loss_orig
-    res_gen_p = 1.0 * gen_p
-    iter_num = 0
-    hydro_constraints = economic_dispatch.make_hydro_constraints_from_res_load_scenario()
-    while True:
-        iter_num += 1
-        load = load_without_loss + all_loss
-        load = pd.DataFrame(load.ravel(), index=datetimes)
-        
-        # "never" decrease (during iteration) some generators
-        min__ = diff_.min()  # this is negative
-        gen_max_pu_t = None
-        gen_min_pu_t = {gen_nm: np.maximum((res_gen_p[:,gen_id] + min__) / economic_dispatch.generators.loc[gen_nm].p_nom,
-                                            env_for_loss.gen_pmin[gen_id] / economic_dispatch.generators.loc[gen_nm].p_nom
-                                            )
-                        for gen_id, gen_nm in enumerate(env_for_loss.name_gen) if env_for_loss.gen_redispatchable[gen_id]}
-        
-        ### run the dispatch with the loss
-        dispatch_res = economic_dispatch.run(load,
-                                             total_solar=total_solar,
-                                             total_wind=total_wind,
-                                             params=params,
-                                             pyomo=False,
-                                             solver_name="cbc",
-                                             gen_constraints=copy.deepcopy(hydro_constraints),
-                                             gen_max_pu_t=gen_max_pu_t,
-                                             gen_min_pu_t=gen_min_pu_t,
-                                             )
-        
-        if dispatch_res is None:     
-            error_ = RuntimeError(f"Pypsa failed to find a solution at iteration {iter_num}")
-            break
-        
-        # assign the generators
-        for gen_id, gen_nm in enumerate(env_for_loss.name_gen):
-            if gen_nm in dispatch_res.chronix.prods_dispatch:
-                res_gen_p[:, gen_id] = 1.0 * dispatch_res.chronix.prods_dispatch[gen_nm].values
-                  
-        #handle wind curtailment
-        mask_winds = total_wind.values > 0.001
-        res_gen_p[mask_winds, :][:,env_for_loss.gen_type == "wind"] *= (dispatch_res.chronix.prods_dispatch['agg_wind'].values[mask_winds] / total_wind.values[mask_winds]).reshape(-1,1)
-        total_wind.loc[mask_winds] = 1.0 * dispatch_res.chronix.prods_dispatch["agg_wind"].values[mask_winds]
-        
-        # re evaluate the losses
-        with warnings.catch_warnings():
-            warnings.filterwarnings("ignore")
-            env_fixed = grid2op.make(
-                env_path,
-                test=True,
-                # grid_path=grid_path, # assign it the 118 grid
-                param=env_param,
-                backend=LightSimBackend(),
-                chronics_class=FromNPY,
-                # chronics_path=path_chronix2grid,
-                data_feeding_kwargs={"load_p": load_p,
-                                     "load_q": load_q,
-                                     "prod_p": 1.0 * res_gen_p,
-                                     "prod_v": gen_v}
-                )
-        diff_ = np.full((env_fixed.max_episode_duration(), env_fixed.n_gen), fill_value=np.NaN)
-        all_loss[:] = np.NaN
-        
-        i = 0
-        obs = env_fixed.reset()
-        all_loss[i] = np.sum(obs.gen_p) - np.sum(obs.load_p)
-        diff_[i] = obs.gen_p - res_gen_p[i]
-
-        
-        done = False
-        while not done:
-            obs, reward, done, info = env_fixed.step(env_fixed.action_space())
-            i += 1
-            if done:
-                break
-            all_loss[i] = np.sum(obs.gen_p) - np.sum(obs.load_p)
-            diff_[i] = obs.gen_p - res_gen_p[i]
-        
-        max_diff_ = np.abs(diff_).max()
-        if not np.isfinite(max_diff_):
-            error_ = RuntimeError(f"Some nans were found in the generated data at iteration {iter_num}")
-            res_gen_p = None
-            quality_ = None
-            break
-            
-            
-        if max_diff_ <= threshold_stop:
-            quality_ = (iter_num,
-                        float(np.mean(np.abs(diff_))),
-                        float(np.percentile(np.abs(diff_), 95)),
-                        float(np.percentile(np.abs(diff_), 99)),
-                        float(max_diff_)
-            )
-            break
-        
-        if iter_num >= iter_quality_decrease:
-            quantile = np.percentile(np.abs(diff_), percentile_quality_decrease)
-            if quantile <= threshold_stop:
-                quality_ = (iter_num,
-                            float(np.mean(np.abs(diff_))),
-                            float(np.percentile(np.abs(diff_), 95)),
-                            float(np.percentile(np.abs(diff_), 99)),
-                            float(np.max(np.abs(diff_)))
-                )
-                break
-                    
-        if iter_num >= max_iter:
-            error_ = RuntimeError("Too much iterations performed when adjusting for the losses")
-            res_gen_p = None
-            quality_ = None
-            break
-        
-    return res_gen_p, error_, quality_
-
-
 def _adjust_gens(all_loss_orig,
                  env_for_loss,
                  datetimes,
@@ -487,6 +254,7 @@ def _adjust_gens(all_loss_orig,
                  max_iter=100,  # declare a failure after this number of iteration
                  iter_quality_decrease=50,  # acept a reduction of the quality after this number of iteration
                  percentile_quality_decrease=99,
+                 hydro_constraints=None,
                  ):
     """This function is an auxilliary function.
     
@@ -542,51 +310,53 @@ def _adjust_gens(all_loss_orig,
     all_loss = all_loss_orig
     res_gen_p = 1.0 * gen_p
     iter_num = 0
-    hydro_constraints = economic_dispatch.make_hydro_constraints_from_res_load_scenario()
+    # hydro_constraints = economic_dispatch.make_hydro_constraints_from_res_load_scenario()
     
     # defined some global variable (used for all optimization problems)
     turned_off_orig = 1.0 * (gen_p[:, env_for_loss.gen_redispatchable] == 0.)
-    ids_hyrdo = []
     total_gen = np.sum(env_for_loss.gen_redispatchable)
     total_step = total_solar.shape[0]
     gen_id = 0
-    for i in range(env_for_loss.n_gen):
-        if env_for_loss.gen_redispatchable[i]:
-            if env_for_loss.gen_type[i] == "hydro":
-                ids_hyrdo.append(gen_id)
-            gen_id += 1
-    ids_hyrdo = np.array(ids_hyrdo)
+    ids_hyrdo = get_gen_ids_hydro(env_for_loss)
     
     # define the constraints    
     scaling_factor = env_for_loss.gen_pmax[env_for_loss.gen_redispatchable]
     p_min = np.repeat(env_for_loss.gen_pmin[env_for_loss.gen_redispatchable].reshape(1,-1)  / scaling_factor,
-                        total_step,
-                        axis=0)
+                      total_step,
+                      axis=0)
     p_max = np.repeat(env_for_loss.gen_pmax[env_for_loss.gen_redispatchable].reshape(1,-1) * params["PmaxErrorCorrRatio"] / scaling_factor,
-                        total_step,
-                        axis=0)
+                      total_step,
+                      axis=0)
+    
+    ramp_min = np.repeat(-env_for_loss.gen_max_ramp_down[env_for_loss.gen_redispatchable].reshape(1,-1) * params["RampErrorCorrRatio"] / scaling_factor,
+                         total_step - 1,
+                         axis=0)
+    ramp_max = np.repeat(env_for_loss.gen_max_ramp_up[env_for_loss.gen_redispatchable].reshape(1,-1) * params["RampErrorCorrRatio"] / scaling_factor,
+                         total_step - 1,
+                         axis=0)
+    
+    if hydro_constraints is not None:
+        p_max[:, ids_hyrdo] = fix_nan_hydro_i_dont_know_why(hydro_constraints["p_max_pu"].values)
+    
+    if "hydro_ramp_reduction_factor" in params:
+        ramp_max[:, ids_hyrdo] /= float(params["hydro_ramp_reduction_factor"])
+        ramp_min[:, ids_hyrdo] /= float(params["hydro_ramp_reduction_factor"])
+        
+    
+    for el, el_nm in zip([p_min, p_max, ramp_min, ramp_max],
+                         ["p_min", "p_max", "ramp_min", "ramp_max"]):
+        if not np.all(np.isfinite(el)):
+            error_ = RuntimeError(f"{el_nm} contains non finite values")
+            return None, error_, None
     
-    ramp_min = np.repeat(-env_for_loss.gen_max_ramp_down[env_for_loss.gen_redispatchable].reshape(1,-1) * params["RampErrorCorrRatio"]  / scaling_factor,
-                            total_step - 1,
-                            axis=0)
-    ramp_max = np.repeat(env_for_loss.gen_max_ramp_up[env_for_loss.gen_redispatchable].reshape(1,-1) * params["RampErrorCorrRatio"]  / scaling_factor,
-                            total_step - 1,
-                            axis=0)
-    p_max[:, ids_hyrdo] = 1.0 * hydro_constraints["p_max_pu"].values
-     
     while True:
         iter_num += 1
-        load = load_without_loss + all_loss
-        load = pd.DataFrame(load.ravel(), index=datetimes)
-        
-        # "never" decrease (during iteration) some generators
-        min__ = diff_.min()  # this is negative
         load = load_without_loss + all_loss - np.sum(res_gen_p[:,~env_for_loss.gen_redispatchable], axis=1)
         scale_for_loads =  np.repeat(scaling_factor.reshape(1,-1), total_step, axis=0)
         target_vector = res_gen_p[:,env_for_loss.gen_redispatchable] / scaling_factor         
         
         #### cvxpy
-        p_t = cp.Variable(shape=(total_step,total_gen), pos=True)
+        p_t = cp.Variable(shape=(total_step,total_gen))
         real_p = cp.multiply(p_t, scale_for_loads)
         
         constraints = [p_t >= p_min,
@@ -598,12 +368,18 @@ def _adjust_gens(all_loss_orig,
         cost = cp.sum_squares(p_t - target_vector) + cp.norm1(cp.multiply(p_t, turned_off_orig))
         prob = cp.Problem(cp.Minimize(cost), constraints)
         try:
-            prob.solve()
+            res = prob.solve()
         except cp.error.SolverError as exc_:
             error_ = RuntimeError(f"cvxpy failed to find a solution at iteration {iter_num}, error {exc_}")
             res_gen_p = None
             quality_ = None
             break
+        
+        if not np.isfinite(res):
+            error_ = RuntimeError(f"cvxpy failed to find a solution at iteration {iter_num}, and return a non finite cost")
+            res_gen_p = None
+            quality_ = None
+            break
             
         # assign the generators
         gen_p_after_optim = real_p.value
@@ -612,7 +388,13 @@ def _adjust_gens(all_loss_orig,
             res_gen_p = None
             quality_ = None
             break
-            
+ 
+        if np.any(~np.isfinite(gen_p_after_optim)):
+            error_ = RuntimeError(f"cvxpy failed to find a solution at iteration {iter_num}, and returned nans.")
+            res_gen_p = None
+            quality_ = None
+            break
+                   
         id_redisp = 0
         for gen_id, gen_nm in enumerate(env_for_loss.name_gen):
             if env_for_loss.gen_redispatchable[gen_id]:
@@ -704,6 +486,7 @@ def _fix_losses_one_scenario(env_for_loss,
                             max_iter=100,  # maximum number of iteration
                             iter_quality_decrease=20,  # after 20 iteration accept a degradation in the quality
                             percentile_quality_decrease=99,  # replace the "at maximum" by "percentile 99%"
+                            hydro_constraints=None,
                             ):
     """This function is an auxilliary function.
     
@@ -743,7 +526,7 @@ def _fix_losses_one_scenario(env_for_loss,
     max_diff_orig = np.zeros(env_for_loss.max_episode_duration())
     datetimes = np.zeros(env_for_loss.max_episode_duration(), dtype=datetime)
     
-    env_for_loss.set_id(scenario_id)
+    # env_for_loss.set_id(scenario_id)  # env_for_loss has only 1 set of time series anyways !
     obs = env_for_loss.reset()
     
     i = 0
@@ -814,7 +597,8 @@ def _fix_losses_one_scenario(env_for_loss,
                                                threshold_stop=threshold_stop,
                                                max_iter=max_iter,
                                                iter_quality_decrease=iter_quality_decrease,
-                                               percentile_quality_decrease=percentile_quality_decrease)
+                                               percentile_quality_decrease=percentile_quality_decrease,
+                                               hydro_constraints=hydro_constraints)
     
     if error_ is not None:
         # the procedure failed
@@ -822,6 +606,32 @@ def _fix_losses_one_scenario(env_for_loss,
     
     return res_gen_p, error_, quality_
 
+def make_env_for_loss(path_env, env_param, load_p, load_q, final_gen_p,
+                      gen_v, start_date_dt, dt_dt):
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore")
+        env_for_loss = grid2op.make(
+            path_env,
+            test=True,
+            # grid_path=grid_path, # assign it the 118 grid
+            param=env_param,
+            backend=LightSimBackend(),
+            chronics_class=FromNPY,
+            # chronics_path=path_chronix2grid,
+            data_feeding_kwargs={"load_p": load_p.values,  # np.concatenate([load_p.values[0].reshape(1,-1), load_p.values]),
+                                    "load_q": load_q.values,  # np.concatenate([load_q.values[0].reshape(1,-1), load_q.values]),
+                                    "prod_p": 1.0 * final_gen_p.values,  # 1.0 * np.concatenate([final_gen_p.values[0].reshape(1,-1), final_gen_p.values]),
+                                    "prod_v": gen_v,  # np.concatenate([gen_v[0].reshape(1,-1), gen_v])}
+                                    "start_datetime": start_date_dt,
+                                    "time_interval": dt_dt,
+            },
+            opponent_budget_per_ts=0.,
+            opponent_init_budget=0.,
+            opponent_class=BaseOpponent,
+            opponent_budget_class=NeverAttackBudget,
+            opponent_action_class=DontAct,
+            )
+    return env_for_loss
 
 def handle_losses(path_env,
                   n_gen,
@@ -839,6 +649,7 @@ def handle_losses(path_env,
                   max_iter=100,
                   iter_quality_decrease=20,  # after 20 iteration accept a degradation in the quality
                   percentile_quality_decrease=99,  # replace the "at maximum" by "percentile 99%"
+                  hydro_constraints=None,
                   ):
     """This function is here to make sure that if you run an AC model with the data generated, then the generator setpoints will not change too much 
     (less than `threshold_stop` MW)
@@ -895,24 +706,9 @@ def handle_losses(path_env,
     env_param.NO_OVERFLOW_DISCONNECTION = True
     gen_v = np.tile(np.array([float(gens_charac.loc[gens_charac["name"] == nm_gen].V) for nm_gen in name_gen ]),
                     load_p.shape[0]).reshape(-1, n_gen)
-    with warnings.catch_warnings():
-        warnings.filterwarnings("ignore")
-        env_for_loss = grid2op.make(
-            path_env,
-            test=True,
-            # grid_path=grid_path, # assign it the 118 grid
-            param=env_param,
-            backend=LightSimBackend(),
-            chronics_class=FromNPY,
-            # chronics_path=path_chronix2grid,
-            data_feeding_kwargs={"load_p": load_p.values,  # np.concatenate([load_p.values[0].reshape(1,-1), load_p.values]),
-                                 "load_q": load_q.values,  # np.concatenate([load_q.values[0].reshape(1,-1), load_q.values]),
-                                 "prod_p": 1.0 * final_gen_p.values,  # 1.0 * np.concatenate([final_gen_p.values[0].reshape(1,-1), final_gen_p.values]),
-                                 "prod_v": gen_v,  # np.concatenate([gen_v[0].reshape(1,-1), gen_v])}
-                                 "start_datetime": start_date_dt,
-                                 "time_interval": dt_dt,
-            }
-            )
+    
+    env_for_loss = make_env_for_loss(path_env, env_param, load_p, load_q,
+                                     final_gen_p, gen_v, start_date_dt, dt_dt)
     res_gen_p, error_, quality_ = _fix_losses_one_scenario(env_for_loss,
                                                            scenario_id,
                                                            loss_param,
@@ -924,15 +720,16 @@ def handle_losses(path_env,
                                                            # after 20 iteration accept a degradation in the quality
                                                            iter_quality_decrease=iter_quality_decrease,  
                                                            # replace the "at maximum" by "percentile 99%"
-                                                           percentile_quality_decrease=percentile_quality_decrease,  
+                                                           percentile_quality_decrease=percentile_quality_decrease,
+                                                           hydro_constraints=hydro_constraints,  
                                                            )
     if error_ is not None:
-        return None, error_, None
+        return None, error_, None, env_for_loss
     
     # reformat the generators
     res_gen_p_df = pd.DataFrame(res_gen_p, index=final_gen_p.index, columns=env_for_loss.name_gen)
     
-    return res_gen_p_df, error_, quality_
+    return res_gen_p_df, error_, quality_, env_for_loss
 
 def save_generated_data(this_scen_path,
                         load_p,
@@ -944,6 +741,7 @@ def save_generated_data(this_scen_path,
                         prod_p,  # after everything, the one the agent sees
                         prod_p_forecasted,
                         debug,
+                        save_load_q=True,
                         sep=';',
                         float_prec=FLOATING_POINT_PRECISION_FORMAT):
     """This function saves the data that have been generated by this script.
@@ -971,13 +769,21 @@ def save_generated_data(this_scen_path,
     float_prec : _type_, optional
         _description_, by default FLOATING_POINT_PRECISION_FORMAT
     """
-    li_dfs = [load_p, load_p_forecasted, load_q, load_q_forecasted, prod_p, prod_p_forecasted]
-    li_nms = ["load_p", "load_p_forecasted", "load_q", "load_q_forecasted", "prod_p", "prod_p_forecasted"]
+    li_dfs = [load_p, load_p_forecasted, prod_p, prod_p_forecasted]
+    li_nms = ["load_p", "load_p_forecasted", "prod_p", "prod_p_forecasted"]
+    
+    if save_load_q:
+        li_dfs.append(load_q)
+        li_nms.append("load_q")
+        li_dfs.append(load_q_forecasted)
+        li_nms.append("load_q_forecasted")
+        
     if debug:
         li_dfs.append(prod_p_generated)
         li_nms.append("prod_p_generated")
         li_dfs.append(prod_p_after_dispatch)
         li_nms.append("prod_p_after_dispatch")
+        
     for df, nm in zip(li_dfs, li_nms):
         df.to_csv(os.path.join(this_scen_path, f'{nm}.csv.bz2'),
                   sep=sep,
@@ -1003,7 +809,14 @@ def save_meta_data(this_scen_path,
                    solar_curtailed_opf,
                    generation_time,
                    saving_time,
-                   files_to_copy=("maintenance_meta.json",)):
+                   forecast_generation_time,
+                   forca_t0_errors,
+                   forca_errors,
+                   amount_curtailed_for,
+                   files_to_copy=("maintenance_meta.json",),
+                   load_ref=None,
+                   solar_ref=None,
+                   wind_ref=None):
     """This function saves the "meta data" required for a succesful grid2op run !
 
     Parameters
@@ -1055,6 +868,7 @@ def save_meta_data(this_scen_path,
                                 "iter_num: number of iteration of the loss algorithm",
                                 "generation_time: total time spent to generate these data (in seconds)",
                                 "saving_time: total time spent to save the generated data (in seconds), this excludes the metadata saving time",
+                                "amount_curtailed_for: sum of all the power that has been curtailed for the forecasts"
                                 ),
                        "total_load": float(total_load),
                        "total_gen": float(total_gen),
@@ -1065,6 +879,10 @@ def save_meta_data(this_scen_path,
                        "solar_curtailed_opf": float(solar_curtailed_opf),
                        "generation_time": float(generation_time),
                        "saving_time": float(saving_time),
+                       "forecast_generation_time": float(forecast_generation_time),
+                       "forecasts_t0_errors": {f"{t0}": f"{error_msg}" for t0, error_msg in zip(forca_t0_errors, forca_errors)},
+                       "forecasts_nb_fail": len(forca_t0_errors),
+                       "amount_curtailed_for": float(np.sum(amount_curtailed_for))
                        },
                   fp=f,
                   sort_keys=True,
@@ -1076,6 +894,13 @@ def save_meta_data(this_scen_path,
             shutil.copy(src=src_path, 
                         dst=os.path.join(this_scen_path, fn_))
     
+    # save the ref curve for load, solar and wind
+    if load_ref is not None:
+        np.save(file=os.path.join(this_scen_path, "load_ref.npy"), arr=load_ref)
+    if solar_ref is not None:
+        np.save(file=os.path.join(this_scen_path, "solar_ref.npy"), arr=solar_ref)
+    if wind_ref is not None:
+        np.save(file=os.path.join(this_scen_path, "wind_ref.npy"), arr=wind_ref)
 
 def generate_a_scenario(path_env,
                         name_gen,
@@ -1093,6 +918,10 @@ def generate_a_scenario(path_env,
                         RampErrorCorrRatio=0.95,
                         threshold_stop=0.05,
                         max_iter=100,
+                        files_to_copy=("maintenance_meta.json",),
+                        save_ref_curve=False,
+                        day_lag=6, # TODO 6 because it's 2050
+                        tol_zero=1e-3,
                         debug=True  # TODO more feature !
                         ):
     """This function generates and save the data for a scenario.
@@ -1148,21 +977,42 @@ def generate_a_scenario(path_env,
     number_of_minutes = int((end_date_dt - start_date_dt).total_seconds() // 60)
     gens_charac = pd.read_csv(os.path.join(path_env, "prods_charac.csv"), sep=",")
     
+    forecast_prng = default_rng(gen_p_forecast_seed)
     # conso generation
-    load_p, load_q, load_p_forecasted, load_q_forecasted = generate_loads(path_env,
-                                                                          load_seed,
-                                                                          start_date_dt,
-                                                                          end_date_dt,
-                                                                          dt,
-                                                                          number_of_minutes,
-                                                                          generic_params,
-                                                                          day_lag=6  # TODO 6 because it's 2050
-                                                                          )
+    tmp_ = generate_loads(path_env,
+                          load_seed,
+                          forecast_prng,
+                          start_date_dt,
+                          end_date_dt,
+                          dt,
+                          number_of_minutes,
+                          generic_params,
+                          day_lag=day_lag 
+                          )
+    (new_forecasts, forecasts_params, load_params, loads_charac,
+     load_p, load_q, load_p_forecasted, load_q_forecasted, load_ref) = tmp_
     
     # renewable energy sources generation
-    res_renew = generate_renewable_energy_sources(path_env,renew_seed, start_date_dt, end_date_dt, dt, number_of_minutes, generic_params, gens_charac)
-    prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted = res_renew
+    res_renew = generate_renewable_energy_sources(path_env,
+                                                  renew_seed,
+                                                  start_date_dt,
+                                                  end_date_dt,
+                                                  dt,
+                                                  number_of_minutes,
+                                                  generic_params,
+                                                  gens_charac,
+                                                  tol_zero=tol_zero)
+    prod_solar, prod_solar_forecasted, prod_wind, prod_wind_forecasted, solar_ref, wind_ref, renew_prng = res_renew
 
+    path_extra_wind = os.path.join(path_env, "wind_extra_params.json")
+    if os.path.exists(path_extra_wind) and os.path.isfile(path_extra_wind):
+        with open(path_extra_wind, "r", encoding="utf-8") as f:
+            extra_winds_params = json.load(f)
+        prod_wind_init = prod_wind
+        prod_wind = apply_maintenance_wind_farm(extra_winds_params, prod_wind_init,
+                                                start_date_dt, end_date_dt, dt,
+                                                renew_prng)
+        
     if prod_solar.isna().any().any():
         error_ = RuntimeError("Nan generated in solar data")
         return error_, None, None, None, None, None, None, None
@@ -1179,43 +1029,78 @@ def generate_a_scenario(path_env,
         final_gen_p[str(el)] = np.NaN
     final_gen_p = final_gen_p[name_gen]
     
+    with open(os.path.join(path_env, "params_opf.json"), "r") as f:
+        opf_params = json.load(f)
+    opf_params["start_date"] = start_date_dt
+    opf_params["end_date"] = end_date_dt
+    opf_params["dt"] = int(dt)
+    opf_params["T"] = number_of_minutes
+    opf_params["planned_std"] = float(generic_params["planned_std"])
+    
     # generate economic dispatch
-    res_disp = generate_economic_dispatch(path_env, start_date_dt, end_date_dt, dt, number_of_minutes, generic_params,
-                                          load_p, prod_solar, prod_wind, name_gen, gen_type, scenario_id, final_gen_p, gens_charac)
-    gen_p_after_dispatch, total_wind_curt_opf, total_solar_curt_opf, error_ = res_disp
+    res_disp = generate_economic_dispatch(path_env, start_date_dt, end_date_dt, dt, number_of_minutes,
+                                          generic_params,
+                                          load_p, prod_solar, prod_wind, name_gen, gen_type, scenario_id,
+                                          final_gen_p, gens_charac, opf_params)
+    gen_p_after_dispatch, total_wind_curt_opf, total_solar_curt_opf, hydro_constraints, error_ = res_disp
     
     if error_ is not None:
         # TODO log that !
         return error_, None, None, None, None, None, None, None
     
     # now try to move the generators so that when I run an AC powerflow, the setpoint of generators does not change "too much"
+    n_gen = len(name_gen)
     if handle_loss:
-        n_gen = len(name_gen)
-        res_gen_p_df, error_, quality_ = handle_losses(path_env,
-                                                       n_gen,
-                                                       name_gen,
-                                                       gens_charac,
-                                                       load_p,
-                                                       load_q,
-                                                       gen_p_after_dispatch,
-                                                       start_date_dt,
-                                                       dt_dt,
-                                                       scenario_id, 
-                                                       PmaxErrorCorrRatio=PmaxErrorCorrRatio,
-                                                       RampErrorCorrRatio=RampErrorCorrRatio,
-                                                       threshold_stop=threshold_stop,
-                                                       max_iter=max_iter)
+        res_gen_p_df, error_, quality_, env_for_loss = handle_losses(path_env,
+                                                                     n_gen,
+                                                                     name_gen,
+                                                                     gens_charac,
+                                                                     load_p,
+                                                                     load_q,
+                                                                     gen_p_after_dispatch,
+                                                                     start_date_dt,
+                                                                     dt_dt,
+                                                                     scenario_id, 
+                                                                     PmaxErrorCorrRatio=PmaxErrorCorrRatio,
+                                                                     RampErrorCorrRatio=RampErrorCorrRatio,
+                                                                     threshold_stop=threshold_stop,
+                                                                     max_iter=max_iter,
+                                                                     hydro_constraints=hydro_constraints)
         if error_ is not None:
             # TODO log that !
             return error_, None, None, None, None, None, None, None
     else:
         res_gen_p_df = 1.0 * gen_p_after_dispatch
         quality_ = (-1, float("Nan"), float("Nan"), float("Nan"), float("Nan"))
+        
+        env_param = Parameters()
+        env_param.NO_OVERFLOW_DISCONNECTION = True
+        gen_v = np.tile(np.array([float(gens_charac.loc[gens_charac["name"] == nm_gen].V) for nm_gen in name_gen ]),
+                        load_p.shape[0]).reshape(-1, n_gen)
+        env_for_loss = make_env_for_loss(path_env, env_param,
+                                         load_p, load_q,
+                                         final_gen_p, gen_v,
+                                         start_date_dt, dt_dt)
     
     prng = default_rng(gen_p_forecast_seed)
-    res_gen_p_forecasted_df = res_gen_p_df * prng.lognormal(mean=0.0, sigma=float(generic_params["planned_std"]), size=res_gen_p_df.shape)
-    res_gen_p_forecasted_df = res_gen_p_forecasted_df.shift(-1)
-    res_gen_p_forecasted_df.iloc[-1] = 1.0 * res_gen_p_forecasted_df.iloc[-2]
+    beg_forca = time.perf_counter()
+    tmp_ = generate_forecasts_gen(new_forecasts,
+                                  prng,
+                                  load_p,
+                                  load_p_forecasted,
+                                  res_gen_p_df,
+                                  float(generic_params["planned_std"]),
+                                  env_for_loss,
+                                  hydro_constraints,
+                                  forecasts_params,
+                                  load_params,
+                                  loads_charac,
+                                  gens_charac,
+                                  path_env,
+                                  opf_params,
+                                  )
+    res_gen_p_forecasted_df_res, amount_curtailed_for, t0_errors, errors = tmp_
+    end_forca = time.perf_counter()
     end_ = time.perf_counter()
     if output_dir is not None:
         beg_save = time.perf_counter()
@@ -1230,7 +1115,7 @@ def generate_a_scenario(path_env,
                             final_gen_p,  # generated, before economic dispatch
                             gen_p_after_dispatch,  # generated, after economic dispatch (and possibly curtailment)
                             res_gen_p_df,
-                            res_gen_p_forecasted_df,
+                            res_gen_p_forecasted_df_res,
                             debug=debug)
         total_load = float(load_p.sum().sum())
         total_gen = float(res_gen_p_df.sum().sum())
@@ -1254,7 +1139,17 @@ def generate_a_scenario(path_env,
                        wind_curtailed_losses=float(wind_curtailed_losses * proper_MWh_unit),
                        solar_curtailed_opf=float(total_solar_curt_opf * proper_MWh_unit),
                        generation_time=end_ - beg_,
-                       saving_time=end_save - beg_save
+                       saving_time=end_save - beg_save,
+                       # forecasts information
+                       forecast_generation_time=end_forca - beg_forca,
+                       forca_t0_errors=t0_errors,
+                       forca_errors=errors,
+                       amount_curtailed_for=amount_curtailed_for,
+                       files_to_copy=files_to_copy,
+                       # ref curve
+                       load_ref=load_ref if save_ref_curve else None,
+                       solar_ref=solar_ref if save_ref_curve else None,
+                       wind_ref=wind_ref if save_ref_curve else None,
                        )
         
-    return error_, quality_, load_p, load_p_forecasted, load_q, load_q_forecasted, res_gen_p_df, res_gen_p_forecasted_df
+    return error_, quality_, load_p, load_p_forecasted, load_q, load_q_forecasted, res_gen_p_df, res_gen_p_forecasted_df_res
diff --git a/chronix2grid/kpi/deterministic/kpis.py b/chronix2grid/kpi/deterministic/kpis.py
index 38ec2db..19e718f 100755
--- a/chronix2grid/kpi/deterministic/kpis.py
+++ b/chronix2grid/kpi/deterministic/kpis.py
@@ -85,12 +85,11 @@ def __init__(self, ref_consumption, syn_consumption, ref_dispatch,
             os.mkdir(os.path.join(self.image_repo, 'thermal_load_kpi'))
             os.mkdir(os.path.join(self.image_repo, 'load_kpi'))
 
-
         # Aggregate chronics variables
-        self.ref_agg_conso = ref_consumption.sum(axis=1)
-        self.syn_agg_conso = syn_consumption.sum(axis=1)
-        self.agg_ref_dispatch = ref_dispatch.sum(axis=1)
-        self.agg_syn_dispatch = synthetic_dispatch.sum(axis=1)
+        self.ref_agg_conso = ref_consumption.drop("Time", axis=1).sum(axis=1)
+        self.syn_agg_conso = syn_consumption.drop("Time", axis=1).sum(axis=1)
+        self.agg_ref_dispatch = ref_dispatch.drop("Time", axis=1).sum(axis=1)
+        self.agg_syn_dispatch = synthetic_dispatch.drop("Time", axis=1).sum(axis=1)
         
         # Read grid characteristics
         if prods_charac is not None:
@@ -147,8 +146,8 @@ def plot_barcharts(self, df_ref, df_syn, save_plots = True, path_name = '', titl
             df_syn = df_syn / df_syn.max()
 
         fig, axes = plt.subplots(1, 2, figsize=(17, 5))
-        sns.barplot(df_ref.index, df_ref, ax=axes[0])
-        sns.barplot(df_syn.index, df_syn, ax=axes[1])
+        sns.barplot(x=df_ref.index, y=df_ref, ax=axes[0])
+        sns.barplot(x=df_syn.index, y=df_syn, ax=axes[1])
         axes[0].set_title('Reference '+title_component, size = 9)
         axes[1].set_title('Synthetic '+title_component, size = 9)
 
@@ -167,14 +166,14 @@ def energy_mix(self, save_plots = True):
         """
 
         # Sum of production per generator type
-        ref_prod_per_gen = self.ref_dispatch.sum(axis = 0)
+        ref_prod_per_gen = self.ref_dispatch.drop("Time", axis=1).sum(axis = 0)
         ref_prod_per_gen = pd.DataFrame({"Prod": ref_prod_per_gen.values, "name":ref_prod_per_gen.index})
         ref_prod_per_gen = ref_prod_per_gen.merge(self.prod_charac[["name","type"]], how = 'left',
                                                   on = 'name')
         ref_prod_per_gen = ref_prod_per_gen.groupby('type').sum()
         ref_prod_per_gen = ref_prod_per_gen.sort_index()
 
-        syn_prod_per_gen = self.syn_dispatch.sum(axis=0)
+        syn_prod_per_gen = self.syn_dispatch.drop("Time", axis=1).sum(axis=0)
         syn_prod_per_gen = pd.DataFrame({"Prod": syn_prod_per_gen.values, "name": syn_prod_per_gen.index})
         syn_prod_per_gen = syn_prod_per_gen.merge(self.prod_charac[["name", "type"]], how='left',
                                                   on='name')
@@ -1177,8 +1176,8 @@ def load_kpi(self, save_plots = True):
         ## Plot results
         # By day
         fig, axes = plt.subplots(1, 2, figsize=(17, 5))
-        sns.barplot(conso_day_ref.index, conso_day_ref, ax=axes[0])
-        sns.barplot(conso_day.index, conso_day, ax=axes[1])
+        sns.barplot(x=conso_day_ref.index, y=conso_day_ref, ax=axes[0])
+        sns.barplot(x=conso_day.index, y=conso_day, ax=axes[1])
         axes[0].set_title('Reference - Normalized load per day of week', size=9)
         axes[1].set_title('Synthetic - Normalized load per day of week', size=9)
 
@@ -1189,8 +1188,8 @@ def load_kpi(self, save_plots = True):
         every_nth = 3
 
         fig, axes = plt.subplots(1, 2, figsize=(17, 5))
-        sns.barplot(conso_week_ref.index, conso_week_ref, ax=axes[0])
-        sns.barplot(conso_week.index, conso_week, ax=axes[1])
+        sns.barplot(x=conso_week_ref.index, y=conso_week_ref, ax=axes[0])
+        sns.barplot(x=conso_week.index, y=conso_week, ax=axes[1])
         axes[0].set_title('Normalized load per week of year', size=9)
         for n, label in enumerate(axes[0].xaxis.get_ticklabels()):
             if n % every_nth != 0:
diff --git a/chronix2grid/kpi/preprocessing/pivot_utils.py b/chronix2grid/kpi/preprocessing/pivot_utils.py
index 2965d25..30b8ebf 100755
--- a/chronix2grid/kpi/preprocessing/pivot_utils.py
+++ b/chronix2grid/kpi/preprocessing/pivot_utils.py
@@ -190,7 +190,12 @@ def eco2mix_to_kpi_regional(kpi_input_folder, timestep, prods_charac, loads_char
         # Time formatting
         df['Space'] = ' '
         df['Time'] = df['Date']+df['Space']+df['Heures']
-        df['Time'] = pd.to_datetime(df['Time'], infer_datetime_format=True)
+        try:
+            # try the default format
+            df["Time"] = pd.to_datetime(df["Time"], format="%d/%m/%Y %H:%M")
+        except ValueError:
+            # otherwise try to infer it
+            df['Time'] = pd.to_datetime(df['Time'], infer_datetime_format=True)
         df.set_index('Time', drop=False, inplace=True)
 
         # Production formatting
@@ -304,7 +309,12 @@ def renewableninja_to_kpi(kpi_input_folder, timestep, loads_charac,
         # Time formatting
         df['Space'] = ' '
         df['Time'] = df['Date'] + df['Space'] + df['Heures']
-        df['Time'] = pd.to_datetime(df['Time'], infer_datetime_format=True)
+        try:
+            # try the default format
+            df["Time"] = pd.to_datetime(df["Time"], format="%d/%m/%Y %H:%M")
+        except ValueError:
+            # otherwise try to infer it
+            df['Time'] = pd.to_datetime(df['Time'], infer_datetime_format=True)
         df.set_index('Time', drop=False, inplace=True)
 
         # Correct some strange values in eco2mix data
diff --git a/requirements.txt b/requirements.txt
index 021824d..e28ed98 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -7,10 +7,10 @@ ipython-genutils>=0.2.0
 ipywidgets>=7.5.1
 json5>=0.9.3
 jupyter>=1.0.0
-matplotlib==3.1.3
+matplotlib>=3.1.3
 numba>=0.48.0
-numpy==1.18.3
-pandas==1.0.5
+numpy>=1.18.3
+pandas>=1.0.5
 pandocfilters>=1.4.2
 pathlib>=1.0.1
 plotly>=4.5.2
@@ -23,7 +23,7 @@ pytz>=2019.3
 #pywinpty==0.5.7
 PyUtilib==5.7.3
 requests>=2.23.0
-seaborn==0.10.0
+seaborn>=0.10.0
 scipy>=1.4.1
 widgetsnbextension>=3.5.1
 
diff --git a/setup.py b/setup.py
index c28643e..55f6c4b 100644
--- a/setup.py
+++ b/setup.py
@@ -66,7 +66,8 @@
                         "widgetsnbextension>=3.5.1",
                         "lightsim2grid",
                         "pypsa",
-                        "cvxpy"
+                        "cvxpy",
+                        "scikit-learn"
                         ],
       extras_require = {
                         "optional": [


Test patch:
diff --git a/tests/integration_tests/test_integration.py b/tests/integration_tests/test_integration.py
index 828edfe..a440cc0 100644
--- a/tests/integration_tests/test_integration.py
+++ b/tests/integration_tests/test_integration.py
@@ -164,7 +164,7 @@ def test_integration_l(self):
         path_ref = self.expected_folder_noloss
         files_to_check=['load_p']
         bool = self.check_frames_equal(path_out, path_ref, files_to_check)
-        self.assertTrue(bool)
+        self.assertTrue(bool, f"error for {files_to_check}")
 
     def test_integration_r(self):
         self.modify_gen_prices(self.mu, self.sigma, self.seed_price_noise_noloss, self.gen_types, self.case_noloss)
@@ -186,7 +186,7 @@ def test_integration_r(self):
         path_ref = self.expected_folder_noloss
         files_to_check=['solar_p','wind_p']
         bool = self.check_frames_equal(path_out, path_ref, files_to_check)
-        self.assertTrue(bool)
+        self.assertTrue(bool, f"error for {files_to_check}")
 
     def test_integration_lrt_nolosscorrection(self):
         self.modify_gen_prices(self.mu, self.sigma, self.seed_price_noise_noloss, self.gen_types, self.case_noloss)
@@ -207,7 +207,7 @@ def test_integration_lrt_nolosscorrection(self):
         path_out = os.path.join(self.generation_output_folder_noloss, "Scenario_0")
         path_ref = self.expected_folder_noloss
         bool = self.check_frames_equal(path_out, path_ref, self.files_tocheck)
-        self.assertTrue(bool)
+        self.assertTrue(bool, f"error for {self.files_tocheck}")
 
     def test_integration_lrt_withlosscorrection(self):
         self.modify_gen_prices(self.mu, self.sigma, self.seed_price_noise_loss, self.gen_types, self.case_loss)
@@ -232,7 +232,7 @@ def test_integration_lrt_withlosscorrection(self):
             ###
             bool = self.check_frames_equal(path_out, path_ref, self.files_tocheck)
             # Check that we obtain the right result dataframe
-            self.assertTrue(bool)
+            self.assertTrue(bool, f"error for {self.files_tocheck}")
             # Check that we have raised a UserWarning for ramp up (one among all warnings that have been raised)
             boolvec_types = [issubclass(w_.category, UserWarning) for w_ in w]
             self.assertTrue(np.any(boolvec_types))
@@ -263,7 +263,7 @@ def test_integration_hydro_reduction(self):
             ###
             bool = self.check_frames_equal(path_out, path_ref, self.files_tocheck)
             # Check that we obtain the right result dataframe
-            self.assertTrue(bool)
+            self.assertTrue(bool, f"error for {self.files_tocheck}")
             # Check that we have raised a UserWarning for ramp up (one among all warnings that have been raised)
             #boolvec_types = [issubclass(w_.category, UserWarning) for w_ in w]
             #self.assertTrue(np.any(boolvec_types))