From a90403dd3be2fe1aea5ad28fa3d9a3b2451941ea Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Tue, 23 Jan 2024 14:53:59 +0100 Subject: [PATCH 1/8] adressing issue #76, adressing bug in data generation for the forecasts --- chronix2grid/__init__.py | 2 +- chronix2grid/grid2op_utils/add_data.py | 45 +++++--- chronix2grid/grid2op_utils/gen_utils.py | 32 ++++-- chronix2grid/grid2op_utils/loads_utils.py | 8 +- chronix2grid/grid2op_utils/utils.py | 6 +- docs/conf.py | 8 +- setup.py | 2 +- tests/integration_tests/test_grid2oputils.py | 111 +++++++++++++++++++ 8 files changed, 179 insertions(+), 35 deletions(-) create mode 100644 tests/integration_tests/test_grid2oputils.py diff --git a/chronix2grid/__init__.py b/chronix2grid/__init__.py index f1e7bde1..ff4ef5a4 100755 --- a/chronix2grid/__init__.py +++ b/chronix2grid/__init__.py @@ -6,4 +6,4 @@ # 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) -___version__ = "1.2.0.post1" +___version__ = "1.2.1" diff --git a/chronix2grid/grid2op_utils/add_data.py b/chronix2grid/grid2op_utils/add_data.py index 1a626df9..03fefc2a 100644 --- a/chronix2grid/grid2op_utils/add_data.py +++ b/chronix2grid/grid2op_utils/add_data.py @@ -19,20 +19,21 @@ 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 + save_ref_curve, day_lag, tol_zero, debug, load_weekly_pattern) = 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) + 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, + load_weekly_pattern=load_weekly_pattern) return res_gen @@ -47,7 +48,8 @@ def add_data(env: grid2op.Environment.Environment, save_ref_curve=False, day_lag=6, # TODO 6 because it's 2050 debug=False, - tol_zero=1e-3 + tol_zero=1e-3, + load_weekly_pattern=None, ): """This function adds some data to already existing scenarios. @@ -68,6 +70,18 @@ def add_data(env: grid2op.Environment.Environment, with_loss: ``bool`` Do you make sure that the generated data will not be modified too much when running with grid2op (default = True). Setting it to False will speed up (by quite a lot) the generation process, but will degrade the data quality. + load_weekly_pattern: pd.DataFrame + The pattern used as a reference to generate the loads. + + It sould be a dataframe with 2 columns: `datetime` and `test`. + + In the first column (`datetime`) you should have time stamp (format "%Y-%m-%d %H:%M:%S" + *eg* `2017-01-07 23:55:00`). The second oui should have a number "approximately one" which + gives the relative ratio of demand for the whole grid at this time stamp. + + We only tested this when the data was given at 5 minutes resolution (two consecutive rows are + distant from 5 minutes) and with the equivalent of 2 years of data. It might work + (or not) in other cases... """ # required parameters @@ -126,7 +140,8 @@ def add_data(env: grid2op.Environment.Environment, save_ref_curve, day_lag, tol_zero, - debug + debug, + load_weekly_pattern )) if nb_core == 1: for args in argss: diff --git a/chronix2grid/grid2op_utils/gen_utils.py b/chronix2grid/grid2op_utils/gen_utils.py index 126cbc04..c5d15ba2 100644 --- a/chronix2grid/grid2op_utils/gen_utils.py +++ b/chronix2grid/grid2op_utils/gen_utils.py @@ -166,8 +166,7 @@ def fix_forecast_ramps(nb_h, 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:] - + 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 @@ -300,7 +299,9 @@ def generate_new_gen_forecasts(prng, 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) + 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: @@ -342,19 +343,30 @@ def generate_forecasts_gen(new_forecasts, 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 + loss_before_ramps_for = 1.01 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() + mask_not_enough_gen = total_gen <= loss_before_ramps_for * total_demand + nb_concerned = (mask_not_enough_gen).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) + rep_factor = np.tile(tmp.reshape(-1,1), nb_concerned).T # how to split 1MW on the controlable generators + est_losses_mw = (loss_before_ramps_for * total_demand - total_gen) # opposite of the loss per step + # we increase the controlable generation when there is not enough generation (loss positive) + # res_gen_p_forecasted_df.loc[mask_not_enough_gen, type(env_for_loss).gen_redispatchable] *= est_losses_mw[mask_not_enough_gen].values.reshape(-1,1) * rep_factor + res_gen_p_forecasted_df.loc[mask_not_enough_gen, type(env_for_loss).gen_redispatchable] += est_losses_mw[mask_not_enough_gen].values.reshape(-1,1) * rep_factor + # the above increase can lead to generator above pmax, when this is the case, + # I cut it + gen_pmax = type(env_for_loss).gen_pmax + for gen_id, is_disp in enumerate(type(env_for_loss).gen_redispatchable): + if not is_disp: + continue + res_gen_p_forecasted_df.iloc[:, gen_id] = np.minimum(res_gen_p_forecasted_df.iloc[:, gen_id].values, + gen_pmax[gen_id]) + # and fix the ramps (for all h) (an optimizer is run t by t) tmp_ = fix_forecast_ramps(nb_h, load_p, load_p_forecasted, diff --git a/chronix2grid/grid2op_utils/loads_utils.py b/chronix2grid/grid2op_utils/loads_utils.py index cff13ae6..3c0dc32e 100644 --- a/chronix2grid/grid2op_utils/loads_utils.py +++ b/chronix2grid/grid2op_utils/loads_utils.py @@ -128,7 +128,8 @@ def generate_loads(path_env, number_of_minutes, generic_params, load_q_from_p_coeff_default=0.7, - day_lag=6): + day_lag=6, + load_weekly_pattern=None): """ This function generates the load for each consumption on a grid @@ -177,7 +178,10 @@ def generate_loads(path_env, 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 load_weekly_pattern is None: + load_weekly_pattern = pd.read_csv(os.path.join(ref_pattern_path, "load_weekly_pattern.csv"), sep=",") + else: + load_weekly_pattern = pd.DataFrame(load_weekly_pattern) if new_forecasts: load_p, load_p_forecasted, load_ref_curve = generate_new_loads(load_seed, diff --git a/chronix2grid/grid2op_utils/utils.py b/chronix2grid/grid2op_utils/utils.py index ea42ca60..68dc29a3 100644 --- a/chronix2grid/grid2op_utils/utils.py +++ b/chronix2grid/grid2op_utils/utils.py @@ -922,7 +922,8 @@ def generate_a_scenario(path_env, save_ref_curve=False, day_lag=6, # TODO 6 because it's 2050 tol_zero=1e-3, - debug=True # TODO more feature ! + debug=True, # TODO more feature ! + load_weekly_pattern=None, ): """This function generates and save the data for a scenario. @@ -987,7 +988,8 @@ def generate_a_scenario(path_env, dt, number_of_minutes, generic_params, - day_lag=day_lag + day_lag=day_lag, + load_weekly_pattern=load_weekly_pattern ) (new_forecasts, forecasts_params, load_params, loads_charac, load_p, load_q, load_p_forecasted, load_q_forecasted, load_ref) = tmp_ diff --git a/docs/conf.py b/docs/conf.py index d43449e3..98b5135a 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -19,12 +19,12 @@ # -- Project information ----------------------------------------------------- project = 'ChroniX2Grid' -copyright = '2020, Antoine Marot, Mario Jothy, Nicolas Megel' -author = 'Antoine Marot, Mario Jothy, Nicolas Megel' +copyright = '2020, Antoine Marot, Mario Jothy, Nicolas Megel, Benjamin Donnot' +author = 'Antoine Marot, Mario Jothy, Nicolas Megel, Benjamin Donnot' # The full version, including alpha/beta/rc tags -release = '1.1.0.post1' -version = '1.1' +release = '1.2.1' +version = '1.2' # -- General configuration --------------------------------------------------- diff --git a/setup.py b/setup.py index 38392286..cd8655dd 100644 --- a/setup.py +++ b/setup.py @@ -15,7 +15,7 @@ setup(name='Chronix2Grid', - version='1.2.0.post1', + version='1.2.1', description='A python package to generate "en-masse" chronics for loads and productions (thermal, renewable)', long_description=long_description, long_description_content_type='text/markdown', diff --git a/tests/integration_tests/test_grid2oputils.py b/tests/integration_tests/test_grid2oputils.py new file mode 100644 index 00000000..7e68d397 --- /dev/null +++ b/tests/integration_tests/test_grid2oputils.py @@ -0,0 +1,111 @@ +# Copyright (c) 2024, 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 copy +import os +import json +import unittest +import warnings +import tempfile +import shutil +import numpy as np +import pandas as pd +import datetime +from packaging import version +import grid2op +from grid2op.Chronics import ChangeNothing + +from chronix2grid.getting_started.example.input.generation.patterns import ref_pattern_path + + +class TestGrid2opUtils(unittest.TestCase): + def setUp(self) -> None: + if version.parse(grid2op.__version__) < version.parse("1.9.8"): + # a fix in grid2Op 1.9.8 : the "loads_charac.csv" was not + # part of the data shipped with the package before + self.skipTest(f"grid2op version too old {grid2op.__version__} < 1.9.8") + return super().setUp() + + def test_not_too_high_value_forecasts(self): + seed = 0 + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + with grid2op.make("l2rpn_idf_2023", test=True) as env: + path = env.get_path_env() + tmp_dir = tempfile.TemporaryDirectory() + new_env_path = os.path.join(tmp_dir.name, "l2rpn_idf_2023") + shutil.copytree(path, new_env_path) + shutil.rmtree(os.path.join(new_env_path, "chronics")) + # keep only the first data (not to generate everything) + with open(os.path.join(new_env_path, "scenario_params.json"), "r") as f: + scenario_params = json.load(f) + scenario_params["all_dates"] = scenario_params["all_dates"][:1] + with open(os.path.join(new_env_path, "scenario_params.json"), "w") as f: + json.dump(fp=f, obj=scenario_params) + env = grid2op.make(new_env_path, + chronics_class=ChangeNothing, + **grid2op.Opponent.get_kwargs_no_opponent()) + env.generate_data(load_weekly_pattern=None, nb_year=1, seed=seed, save_ref_curve=True) + gen_p_for_orig = pd.read_csv(os.path.join(new_env_path, + "chronics", + "2035-01-01_0", + "prod_p_forecasted.csv.bz2"), + sep=";") + assert (gen_p_for_orig.iloc[:,2] <= type(env).gen_pmax[2]).all() + tmp_dir.cleanup() + + def test_load_weekly_pattern(self): + seed = 0 + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + with grid2op.make("l2rpn_wcci_2022", test=True) as env: + path = env.get_path_env() + tmp_dir = tempfile.TemporaryDirectory() + new_env_path = os.path.join(tmp_dir.name, "l2rpn_wcci_2022") + shutil.copytree(path, new_env_path) + shutil.rmtree(os.path.join(new_env_path, "chronics")) + # keep only the first data (not to generate everything) + with open(os.path.join(new_env_path, "scenario_params.json"), "r") as f: + scenario_params = json.load(f) + scenario_params["all_dates"] = scenario_params["all_dates"][:1] + with open(os.path.join(new_env_path, "scenario_params.json"), "w") as f: + json.dump(fp=f, obj=scenario_params) + env = grid2op.make(new_env_path, + chronics_class=ChangeNothing, + **grid2op.Opponent.get_kwargs_no_opponent()) + env.generate_data(load_weekly_pattern=None, nb_year=1, seed=seed, save_ref_curve=True) + load_weekly_pattern = pd.read_csv(os.path.join(ref_pattern_path, "load_weekly_pattern.csv"), sep=",") + load_ref_orig = np.load(os.path.join(new_env_path, "chronics", "2050-01-03_0", "load_ref.npy")) + total_demand_orig = np.sum(load_ref_orig, axis=1) + + # change the load weekly pattern + load_weekly_pattern2 = load_weekly_pattern[load_weekly_pattern["datetime"] >= "2018"] + load_weekly_pattern3 = load_weekly_pattern[load_weekly_pattern["datetime"] < "2018"] + load_weekly_pattern_new = pd.concat([load_weekly_pattern2, load_weekly_pattern3]) + load_weekly_pattern_new.reset_index(inplace=True) + load_weekly_pattern_new = load_weekly_pattern_new[["datetime", "test"]] + load_weekly_pattern_new["datetime"] = copy.deepcopy(load_weekly_pattern["datetime"]) + # delete original data + shutil.rmtree(os.path.join(new_env_path, "chronics")) + + # start a new generation + env.generate_data(load_weekly_pattern=load_weekly_pattern_new, nb_year=1, seed=seed, save_ref_curve=True) + load_ref_new = np.load(os.path.join(new_env_path, "chronics", "2050-01-03_0", "load_ref.npy")) + total_demand_new = np.sum(load_ref_new, axis=1) + + # recompute ref case to make sure it works + shutil.rmtree(os.path.join(new_env_path, "chronics")) + env.generate_data(load_weekly_pattern=None, nb_year=1, seed=seed, save_ref_curve=True) + load_ref_orig2 = np.load(os.path.join(new_env_path, "chronics", "2050-01-03_0", "load_ref.npy")) + total_demand_orig2 = np.sum(load_ref_orig2, axis=1) + + # compate the ref curves and the load data + assert np.allclose(total_demand_orig, total_demand_orig2) + assert not np.allclose(total_demand_orig, total_demand_new) + tmp_dir.cleanup() + \ No newline at end of file From 8c8484fd324d6c8b4bd618d66760dfd6ed0ff0e1 Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Wed, 3 Apr 2024 16:54:28 +0200 Subject: [PATCH 2/8] fixing some issue, adding some warnings and exceptions when errors are detected, improved the API with grid2op --- .../PypsaEconomicDispatch.py | 258 ++++++++++++------ .../run_economic_dispatch.py | 47 +++- .../_EDispatch_L2RPN2020/utils.py | 51 +++- .../consumption/consumption_utils.py | 14 +- .../renewable/generate_solar_wind.py | 6 +- .../generation/renewable/solar_wind_utils.py | 14 +- chronix2grid/grid2op_utils/add_data.py | 38 ++- chronix2grid/grid2op_utils/utils.py | 166 ++++++++--- 8 files changed, 438 insertions(+), 156 deletions(-) diff --git a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/PypsaEconomicDispatch.py b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/PypsaEconomicDispatch.py index c4ed0815..c64f4a94 100644 --- a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/PypsaEconomicDispatch.py +++ b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/PypsaEconomicDispatch.py @@ -38,15 +38,28 @@ class PypsaDispatcher(Dispatcher, pypsa.Network): # PATCH # to avoid problems for respecting pmax and ramps when rounding production values in chronics at the end, we apply a correcting factor - PmaxCorrectingFactor = 1 + PmaxCorrectingFactor = 1. RampCorrectingFactor = 0.1 - - def __init__(self, *args, **kwargs): + DebugPmax = 1e9 # in debug mode, pmax for the "infinite" generator + DebugCost = 1e7 # in debug mode, cost for the "infinite" generator + + def __init__(self, + *args, + pypsa_debug_add_inf_gen=False, + pypsa_debug_add_thermal=True, + pypsa_debug_add_hydro=True, + pypsa_debug_add_nuc=True, + pypsa_debug_add_solar=True, + pypsa_debug_add_wind=True, + pypsa_debug_add_slack=True, + **kwargs): super().__init__(*args, **kwargs) - self.add('Bus', 'node') - self.add('Load', name='agg_load', bus='node') + self.add('Bus', 'only_bus') + self.add('Load', name='agg_load', bus='only_bus') self._env = None # The grid2op environment when instanciated with from_gri2dop_env self._df = None + + self._all_gen_names = set() self._chronix_scenario = None self._simplified_chronix_scenario = None self._has_results = False @@ -54,9 +67,103 @@ def __init__(self, *args, **kwargs): self._pmax_solar = None self._pmax_wind = None + self._pypsa_debug_add_inf_gen = bool(pypsa_debug_add_inf_gen) + if self._pypsa_debug_add_inf_gen: + self.add(class_name='Generator', + name="debug_infinite_gen", + bus='only_bus', + marginal_cost=type(self).DebugCost, + p_nom=type(self).DebugPmax, + ramp_limit_up=1., + ramp_limit_down=1., + ) + self._pypsa_debug_add_thermal = bool(pypsa_debug_add_thermal) + self._pypsa_debug_add_hydro = bool(pypsa_debug_add_hydro) + self._pypsa_debug_add_nuc = bool(pypsa_debug_add_nuc) + self._pypsa_debug_add_solar = bool(pypsa_debug_add_solar) + self._pypsa_debug_add_wind = bool(pypsa_debug_add_wind) + self._pypsa_debug_add_slack = bool(pypsa_debug_add_slack) + self._pypsa_debug_flags = {"wind": self._pypsa_debug_add_wind, + "solar": self._pypsa_debug_add_solar, + "hydro": self._pypsa_debug_add_hydro, + "thermal": self._pypsa_debug_add_thermal, + "nuclear": self._pypsa_debug_add_nuc, + "slack": self._pypsa_debug_add_slack, + "infinite_gen": self._pypsa_debug_add_inf_gen, + "slack_maybe_ignored": ((not self._pypsa_debug_add_nuc) or + (not self._pypsa_debug_add_thermal) or + (not self._pypsa_debug_add_hydro) or + (not self._pypsa_debug_add_slack))} @classmethod - def from_gri2op_env(cls, grid2op_env): + def _aux_add_solar(cls, net): + if net._pypsa_debug_add_solar and net._pmax_solar > 0.: + net.add('Generator', + name='agg_solar', + bus='only_bus', + carrier="solar", + marginal_cost=0., + p_nom=net._pmax_solar + ) + elif not net._pypsa_debug_add_solar and net._pmax_solar > 0.: + print(f"Solar generators are ignored in the OPF due to a pysa_debug flag") + + @classmethod + def _aux_add_wind(cls, net): + if net._pypsa_debug_add_wind and net._pmax_wind > 0.: + net.add('Generator', + name='agg_wind', + bus='only_bus', + carrier="wind", + marginal_cost=0.1, # we prefer to curtail the wind if we have the choice + # that's because solar should be distributed on the grid + p_nom=net._pmax_wind + ) + elif not net._pypsa_debug_add_wind and net._pmax_wind > 0.: + print(f"Solar generators are ignored in the OPF due to a pysa_debug flag") + + @classmethod + def _aux_add_controlable_gen(cls, net, gen_type, carrier_types_to_exclude, p_max, name, ramp_up, ramp_down, gen_cost_per_MW): + net._all_gen_names.add(name) + if gen_type not in carrier_types_to_exclude: + pnom = p_max - cls.PmaxCorrectingFactor + if pnom <= 0.: + print(f"\t Issue for {name}: a negative pmax would be applied. We put {cls.PmaxCorrectingFactor} instead") + pnom = max(pnom, cls.PmaxCorrectingFactor) + + rampUp = (ramp_up - cls.RampCorrectingFactor) / p_max + rampDown = (ramp_down - cls.RampCorrectingFactor) / p_max + if rampUp <= 0.: + print(f"\t Issue for {name}: a negative ramp_up would be applied. We put {cls.RampCorrectingFactor / p_max} instead") + rampUp = cls.RampCorrectingFactor / p_max + if rampDown <= 0.: + print(f"\t Issue for {name}: a negative ramp_down would be applied. We put {cls.RampCorrectingFactor / p_max} instead") + rampDown = cls.RampCorrectingFactor / p_max + if net._pypsa_debug_flags[gen_type]: + net.add( + class_name='Generator', + name=name, + bus='only_bus', + p_nom=pnom, + carrier=gen_type, + marginal_cost=gen_cost_per_MW, + ramp_limit_up=rampUp, + ramp_limit_down=rampDown, + ) + else: + print(f"\t generator {name} (type {gen_type}) is ignored because of a pypsa_debug flag is set for its type") + + @classmethod + def from_gri2op_env(cls, + grid2op_env, + *, + pypsa_debug_add_inf_gen=False, + pypsa_debug_add_thermal=True, + pypsa_debug_add_hydro=True, + pypsa_debug_add_nuc=True, + pypsa_debug_add_solar=True, + pypsa_debug_add_wind=True, + pypsa_debug_add_slack=True,): """ Implements the abstract method of *Dispatcher* @@ -68,54 +175,49 @@ def from_gri2op_env(cls, grid2op_env): ------- net: :class:`pypsa.Network` """ - net = cls() + net = cls(pypsa_debug_add_inf_gen=pypsa_debug_add_inf_gen, + pypsa_debug_add_thermal=pypsa_debug_add_thermal, + pypsa_debug_add_hydro=pypsa_debug_add_hydro, + pypsa_debug_add_nuc=pypsa_debug_add_nuc, + pypsa_debug_add_solar=pypsa_debug_add_solar, + pypsa_debug_add_wind=pypsa_debug_add_wind, + pypsa_debug_add_slack=pypsa_debug_add_slack, + ) net._env = grid2op_env - carrier_types_to_exclude = ['wind', 'solar'] - for i, generator in enumerate(grid2op_env.name_gen): - gen_type = grid2op_env.gen_type[i] - if gen_type not in carrier_types_to_exclude: - p_max = grid2op_env.gen_pmax[i] - pnom = p_max - cls.PmaxCorrectingFactor - - rampUp = (grid2op_env.gen_max_ramp_up[i] - cls.RampCorrectingFactor) / p_max - rampDown = (grid2op_env.gen_max_ramp_down[i] - cls.RampCorrectingFactor) / p_max - - net.add( - class_name='Generator', name=generator, bus='node', - p_nom=pnom, carrier=grid2op_env.gen_type[i], - marginal_cost=grid2op_env.gen_cost_per_MW[i], - ramp_limit_up=rampUp, - ramp_limit_down=rampDown, - ) - # add total wind and solar (for curtailment) - net._pmax_solar = np.sum([grid2op_env.gen_pmax[i] - for i in range(grid2op_env.n_gen) - if grid2op_env.gen_type[i] == "solar"]) - net.add('Generator', - name='agg_solar', - bus='node', - carrier="solar", - marginal_cost=0., - p_nom=net._pmax_solar - ) + net._pmax_solar = np.sum([env_cls.gen_pmax[i] + for i in range(env_cls.n_gen) + if env_cls.gen_type[i] == "solar"]) + cls._aux_add_solar(net) + net._pmax_wind = np.sum([grid2op_env.gen_pmax[i] - for i in range(grid2op_env.n_gen) - if grid2op_env.gen_type[i] == "wind"]) - net.add('Generator', - name='agg_wind', - bus='node', - carrier="wind", - marginal_cost=0.1, # we prefer to curtail the wind if we have the choice - # that's because solar should be distributed on the grid - p_nom=net._pmax_wind - ) + for i in range(grid2op_env.n_gen) + if grid2op_env.gen_type[i] == "wind"]) + cls._aux_add_wind(net) + carrier_types_to_exclude = ['wind', 'solar'] + env_cls = type(grid2op_env) + for i, generator in enumerate(env_cls.name_gen): + gen_type = env_cls.gen_type[i] + p_max = env_cls.gen_pmax[i] + ramp_up = env_cls.gen_max_ramp_up[i] + ramp_down = env_cls.gen_max_ramp_down[i] + gen_cost_per_MW = env_cls.gen_cost_per_MW[i] + cls._aux_add_controlable_gen(net, gen_type, carrier_types_to_exclude, p_max, generator, ramp_up, ramp_down, gen_cost_per_MW) return net @classmethod - def from_dataframe(cls, env_df): + def from_dataframe(cls, + env_df, + *, + pypsa_debug_add_inf_gen=False, + pypsa_debug_add_thermal=True, + pypsa_debug_add_hydro=True, + pypsa_debug_add_nuc=True, + pypsa_debug_add_solar=True, + pypsa_debug_add_wind=True, + pypsa_debug_add_slack=True): """ Implements the abstract method of *Dispatcher* @@ -130,60 +232,39 @@ def from_dataframe(cls, env_df): with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=FutureWarning) # issue with pypsa on some version of numpy - net = cls() + net = cls(pypsa_debug_add_inf_gen=pypsa_debug_add_inf_gen, + pypsa_debug_add_thermal=pypsa_debug_add_thermal, + pypsa_debug_add_hydro=pypsa_debug_add_hydro, + pypsa_debug_add_nuc=pypsa_debug_add_nuc, + pypsa_debug_add_solar=pypsa_debug_add_solar, + pypsa_debug_add_wind=pypsa_debug_add_wind, + pypsa_debug_add_slack=pypsa_debug_add_slack, + ) net._df = env_df - - carrier_types_to_exclude = ['wind', 'solar'] - - for i, (generator, gen_type, p_max, ramp_up, ramp_down, gen_cost_per_MW) in enumerate(zip(env_df['name'], - env_df['type'], - env_df['pmax'], - env_df['max_ramp_up'], - env_df['max_ramp_down'], - env_df['cost_per_mw'])): - if gen_type not in carrier_types_to_exclude: - pnom = (p_max - cls.PmaxCorrectingFactor) - rampUp = (ramp_up- cls.RampCorrectingFactor) / p_max - rampDown = (ramp_down - cls.RampCorrectingFactor) / p_max - - net.add( - class_name='Generator', - name=generator, - bus='node', - p_nom=pnom, - carrier=gen_type, - marginal_cost=gen_cost_per_MW, - ramp_limit_up=rampUp, - ramp_limit_down=rampDown, - ) - + # add total wind and solar (for curtailment) net._pmax_solar = np.sum([p_max for i, (gen_type, p_max) in enumerate(zip(env_df['type'], env_df['pmax'])) if gen_type == "solar"]) - net.add('Generator', - name='agg_solar', - bus='node', - carrier="solar", - marginal_cost=0., - p_nom=net._pmax_solar - ) + cls._aux_add_solar(net) net._pmax_wind = np.sum([p_max for i, (gen_type, p_max) in enumerate(zip(env_df['type'], env_df['pmax'])) if gen_type == "wind"]) - net.add('Generator', - name='agg_wind', - bus='node', - carrier="wind", - marginal_cost=0.1, # we prefer to curtail the wind if we have the choice - # that's because solar should be distributed on the grid - p_nom=net._pmax_wind - ) - return net + cls._aux_add_wind(net) + + carrier_types_to_exclude = ['wind', 'solar'] + for i, (generator, gen_type, p_max, ramp_up, ramp_down, gen_cost_per_MW) in enumerate(zip(env_df['name'], + env_df['type'], + env_df['pmax'], + env_df['max_ramp_up'], + env_df['max_ramp_down'], + env_df['cost_per_mw'])): + cls._aux_add_controlable_gen(net, gen_type, carrier_types_to_exclude, p_max, generator, ramp_up, ramp_down, gen_cost_per_MW) + return net def run(self, load, @@ -207,6 +288,7 @@ def run(self, total_solar = total_solar / self._pmax_solar if total_wind is not None: total_wind = total_wind / self._pmax_wind + prods_dispatch, terminal_conditions, marginal_prices = \ main_run_disptach( self if not by_carrier else self.simplify_net(), diff --git a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/run_economic_dispatch.py b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/run_economic_dispatch.py index 252f7e38..455dc4da 100644 --- a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/run_economic_dispatch.py +++ b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/run_economic_dispatch.py @@ -9,6 +9,7 @@ import argparse import os import time +from typing import Dict, Literal import pandas as pd import pypsa @@ -37,8 +38,11 @@ def main_run_disptach(pypsa_net, # values passed by the users and params if gen_constraints is None: gen_constraints = {} - gen_constraints = update_gen_constrains(gen_constraints) - + + gen_constraints : Dict[Literal["p_min_pu", "p_max_pu"], pd.DataFrame]= update_gen_constrains(gen_constraints) + for el, df in gen_constraints.items(): + df.bfill(inplace=True) # propagate possible NaNs "backward", good for first Na values + df.ffill(inplace=True) # propagate possible NaNs "forward", good if last Na values params = update_params(load.shape[0], load.index[0], params) print ('Preprocessing input data..') @@ -73,16 +77,43 @@ def main_run_disptach(pypsa_net, slack_name = None slack_pmin = None slack_pmax = None - if 'slack_name' in params: - if 'slack_pmin' in params: + slack_ignored_debug_flag = False # is the slack ignored due to a debug flag + if 'slack_name' in params and pypsa_net._pypsa_debug_add_slack: + slack_name = str(params["slack_name"]) + if slack_name not in pypsa_net.generators.index: + if pypsa_net._pypsa_debug_flags["slack_maybe_ignored"]: + if slack_name not in pypsa_net._all_gen_names: + raise RuntimeError("Wrong configuration detected: make sure that the slack is a generator " + "but also a controlable one !") + slack_ignored_debug_flag = True # in this case it is ! + else: + print(pypsa_net.generators.index) + raise RuntimeError("Wrong configuration detected: make sure that the slack is a generator " + "but also a controlable one !") + + if 'slack_pmin' in params and not slack_ignored_debug_flag: # user specified a pmin for the slack bus - slack_name = str(params["slack_name"]) slack_pmin = float(params["slack_pmin"]) / float(pypsa_net.generators.loc[slack_name].p_nom) - - if 'slack_pmax' in params: + if slack_name not in gen_constraints["p_min_pu"]: + print(gen_constraints["p_min_pu"].columns) + raise RuntimeError("Unknown error in chronix2grid: the provided slack has no constraints on the pypsa grid...") + if slack_pmin < gen_constraints["p_min_pu"][slack_name].min(): + raise RuntimeError("Wrong configuration detected: you put in 'params_opf' " + "a minimum value for the slack ('slack_pmin') below its pmin.") + if slack_pmin > 1.: + raise RuntimeError("Unfeasible parameters: slack pmin would be higher than its pmax.") + if 'slack_pmax' in params and not slack_ignored_debug_flag: # user specified a pmin for the slack bus - slack_name = str(params["slack_name"]) slack_pmax = float(params["slack_pmax"]) / float(pypsa_net.generators.loc[slack_name].p_nom) + if slack_name not in gen_constraints["p_max_pu"]: + print(gen_constraints["p_min_pu"].columns) + raise RuntimeError("Unknown error in chronix2grid: the provided slack has no constraints on the pypsa grid...") + if slack_pmax > gen_constraints["p_max_pu"][slack_name].max(): + raise RuntimeError("Wrong configuration detected: you put in 'params_opf' " + "a maximum value for the slack ('slack_pmax') above its pmax.") + if slack_pmax > 1.: + print("Doubtful parameters: slack pmax would be higher than its pmax, basically " + "the parameters `slack_pmin` of `params_opf.json` is ignored.") error_ = False start = time.time() diff --git a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py index c0f042db..cfb209aa 100644 --- a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py +++ b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py @@ -389,8 +389,18 @@ def run_opf(net, gen_max[slack_name] = slack_pmax if slack_name is not None and "slack_ramp_limit_ratio" in params: - net.generators.ramp_limit_up[slack_name] *= float(params["slack_ramp_limit_ratio"]) - net.generators.ramp_limit_down[slack_name] *= float(params["slack_ramp_limit_ratio"]) + if slack_name in net.generators.ramp_limit_up: + net.generators.ramp_limit_up[slack_name] *= float(params["slack_ramp_limit_ratio"]) + else: + if not net._pypsa_debug_flags["slack_maybe_ignored"]: + # slack should be there ! + raise RuntimeError("Slack is missing from the grid generator but should be there") + if slack_name in net.generators.ramp_limit_down: + net.generators.ramp_limit_down[slack_name] *= float(params["slack_ramp_limit_ratio"]) + else: + if not net._pypsa_debug_flags["slack_maybe_ignored"]: + # slack should be there ! + raise RuntimeError("Slack is missing from the grid generator but should be there") if gen_max_pu_t is not None: # addition contraint on the max_pu, used for example when splitting the loss @@ -412,6 +422,43 @@ def run_opf(net, net.generators_t.p_max_pu = pd.concat([gen_max], axis=1) net.generators_t.p_min_pu = pd.concat([gen_min], axis=1) + # raise error for non feasible configurations + if (net.generators["p_min_pu"].values < 0.).any(): + raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative p_min for {net.generators.loc[net.generators['p_min_pu'].values < 0.].index}") + if (net.generators["p_max_pu"].values < 0.).any(): + raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative p_max for {net.generators.loc[net.generators['p_max_pu'].values < 0.].index}") + if (net.generators["ramp_limit_up"].values < 0.).any(): + raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative ramps (up) for {net.generators.loc[net.generators['ramp_limit_up'].values < 0.].index}") + if (net.generators["ramp_limit_down"].values < 0.).any(): + raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative ramps (down) for {net.generators.loc[net.generators['ramp_limit_down'].values < 0.].index}") + if (net.generators_t.p_min_pu.min() > 1.).any(): + debug_info = net.generators_t.p_min_pu.min() + raise RuntimeError(f"A generator would have a p_min above its pmax ! most likely the slack... Check {debug_info.loc[debug_info > 1.].index.values}") + if ((net.generators_t["p_max_pu"] - net.generators_t["p_min_pu"]) < 0.).any().any(): + tmp_debug = ((net.generators_t["p_max_pu"] - net.generators_t["p_min_pu"]) < 0.).any() + raise RuntimeError(f"Invalid configuration: infeasibility detected for some generators (most likely the slack), check generator {tmp_debug[tmp_debug].index.values}") + + # try to detect some non feasible states + if "agg_wind" in net.generators_t["p_max_pu"]: + total_wind = (net.generators.loc["agg_wind", "p_nom"] * net.generators_t["p_max_pu"]["agg_wind"]).values + else: + total_wind = 0. + if "agg_solar" in net.generators_t["p_max_pu"]: + total_solar = (net.generators.loc["agg_solar", "p_nom"] * net.generators_t["p_max_pu"]["agg_solar"]).values + else: + total_solar = 0. + total_pmax = net.generators.loc[(net.generators.index != "agg_wind") & (net.generators.index != "agg_solar"), "p_nom"].sum() + tmp_ = net.generators.loc[(net.generators.index != "agg_wind") & (net.generators.index != "agg_solar"), ["p_nom", "p_min_pu"]].values + total_pmin = (tmp_[:, 0] * tmp_[:, 1]).sum() + min_net_demand = demand["agg_load"].values - total_solar - total_solar + if (min_net_demand > total_pmax).any(): + return None, RuntimeError("Some non feasible time steps are found (net demand above pmax of controlable generators)") + if (demand["agg_load"].values < total_pmin).any(): + return None, RuntimeError("Some non feasible time steps are found (demand [with full curtailment of solar and wind] is below pmin of controlable generators)") + + # TODO check also feasible ramps + # tmp_ = net.generators.loc[(net.generators.index != "agg_wind") & (net.generators.index != "agg_solar"), ["p_nom", "ramp_limit_up", "ramp_limit_down"]].values + # ++ ++ ++ ++ # Run Linear OPF status, termination_condition = net.lopf(net.snapshots, **kwargs) diff --git a/chronix2grid/generation/consumption/consumption_utils.py b/chronix2grid/generation/consumption/consumption_utils.py index 9f772e85..2acbc504 100755 --- a/chronix2grid/generation/consumption/consumption_utils.py +++ b/chronix2grid/generation/consumption/consumption_utils.py @@ -182,10 +182,16 @@ def create_csv(prng, dict_, path, forecasted=False, reordering=True, noise=None, df = df.head(len(df)-1) # Last value is lonely for another day if reordering: value = [] - for name in list(df): - value.append(utils.natural_keys(name)) - new_ordering = [x for _, x in sorted(zip(value, list(df)))] - df = df[new_ordering] + grid2op_default_name = True + try: + for name in list(df): + value.append(utils.natural_keys(name)) + except IndexError: + grid2op_default_name = False + # cannot reorder in this case + if grid2op_default_name: + new_ordering = [x for _ ,x in sorted(zip(value ,list(df)))] + df = df[new_ordering] if shift: df = df.shift(-1) df = df.fillna(0) diff --git a/chronix2grid/generation/renewable/generate_solar_wind.py b/chronix2grid/generation/renewable/generate_solar_wind.py index d4a18d3b..3e06f1bb 100755 --- a/chronix2grid/generation/renewable/generate_solar_wind.py +++ b/chronix2grid/generation/renewable/generate_solar_wind.py @@ -204,7 +204,8 @@ def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, prod_wind = swutils.create_csv( prng, - wind_series, os.path.join(scenario_destination_path, 'wind_p.csv.bz2') if scenario_destination_path is not None else None, + wind_series, + os.path.join(scenario_destination_path, 'wind_p.csv.bz2') if scenario_destination_path is not None else None, reordering=True, noise=params['planned_std'], write_results=write_results @@ -212,7 +213,8 @@ def main(scenario_destination_path, seed, params, prods_charac, solar_pattern, prod_p = swutils.create_csv( prng, - prods_series, os.path.join(scenario_destination_path, 'prod_p.csv.bz2') if scenario_destination_path is not None else None, + prods_series, + os.path.join(scenario_destination_path, 'prod_p.csv.bz2') if scenario_destination_path is not None else None, reordering=True, noise=params['planned_std'], write_results=write_results diff --git a/chronix2grid/generation/renewable/solar_wind_utils.py b/chronix2grid/generation/renewable/solar_wind_utils.py index 60e8b611..8ea7a137 100755 --- a/chronix2grid/generation/renewable/solar_wind_utils.py +++ b/chronix2grid/generation/renewable/solar_wind_utils.py @@ -242,10 +242,16 @@ def create_csv(prng, dict_, path, reordering=True, noise=None, shift=False, df = df.head(len(df ) -1) if reordering: value = [] - for name in list(df): - value.append(utils.natural_keys(name)) - new_ordering = [x for _ ,x in sorted(zip(value ,list(df)))] - df = df[new_ordering] + grid2op_default_name = True + try: + for name in list(df): + value.append(utils.natural_keys(name)) + except IndexError: + grid2op_default_name = False + # cannot reorder in this case + if grid2op_default_name: + new_ordering = [x for _ ,x in sorted(zip(value ,list(df)))] + df = df[new_ordering] if noise is not None: df *= ( 1 +noise * prng.normal(0, 1, df.shape)) #df *= (1 + noise * np.random.normal(0, 1, df.shape)) #older version - to be removed diff --git a/chronix2grid/grid2op_utils/add_data.py b/chronix2grid/grid2op_utils/add_data.py index 03fefc2a..2b434541 100644 --- a/chronix2grid/grid2op_utils/add_data.py +++ b/chronix2grid/grid2op_utils/add_data.py @@ -13,27 +13,32 @@ from chronix2grid.grid2op_utils.utils import generate_a_scenario, get_last_scenario_id from numpy.random import default_rng - + # 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, load_weekly_pattern) = 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, + save_ref_curve, day_lag, tol_zero, debug, load_weekly_pattern, + kwargs) = args + res_gen = generate_a_scenario(path_env=path_env, + name_gen=name_gen, + gen_type=gen_type, + output_dir=output_dir, + start_date=start_date, + dt=dt, + scen_id=scen_id, + load_seed=load_seed, + renew_seed=renew_seed, + gen_p_forecast_seed=gen_p_forecast_seed, + handle_loss=handle_loss, files_to_copy=files_to_copy, save_ref_curve=save_ref_curve, day_lag=day_lag, tol_zero=tol_zero, debug=debug, - load_weekly_pattern=load_weekly_pattern) + load_weekly_pattern=load_weekly_pattern, + **kwargs) return res_gen @@ -50,6 +55,7 @@ def add_data(env: grid2op.Environment.Environment, debug=False, tol_zero=1e-3, load_weekly_pattern=None, + **kwargs ): """This function adds some data to already existing scenarios. @@ -118,7 +124,6 @@ def add_data(env: grid2op.Environment.Environment, path_env = env.get_path_env() name_gen = env.name_gen gen_type = env.gen_type - errors = {} argss = [] for j, scen_id in enumerate(scen_ids): for i, start_date in enumerate(li_months): @@ -141,18 +146,21 @@ def add_data(env: grid2op.Environment.Environment, day_lag, tol_zero, debug, - load_weekly_pattern + load_weekly_pattern, + kwargs, )) if nb_core == 1: + errors = {} for args in argss: res_gen = generate_a_scenario_wrapper(args) error_, *_ = res_gen if error_ is not None: + error_date = args[4] print("=============================") - print(f" Error for {start_date} {scen_id} ") + print(f" Error for {error_date} {scen_id} ") print(f"{error_}") print("=============================") - errors[f'{start_date}_{scen_id}'] = f"{error_}" + errors[f'{error_date}_{scen_id}'] = f"{error_}" # load previous data path_json_error = os.path.join(output_dir, "errors.json") diff --git a/chronix2grid/grid2op_utils/utils.py b/chronix2grid/grid2op_utils/utils.py index 68dc29a3..07b7885e 100644 --- a/chronix2grid/grid2op_utils/utils.py +++ b/chronix2grid/grid2op_utils/utils.py @@ -139,7 +139,15 @@ def generate_renewable_energy_sources(path_env, renew_seed, start_date_dt, end_d 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, - opf_params): + opf_params, + pypsa_debug_add_inf_gen=False, + pypsa_debug_add_thermal=True, + pypsa_debug_add_hydro=True, + pypsa_debug_add_nuc=True, + pypsa_debug_add_solar=True, + pypsa_debug_add_wind=True, + pypsa_debug_add_slack=True, + ): """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 @@ -179,6 +187,7 @@ def generate_economic_dispatch(path_env, start_date_dt, end_date_dt, dt, number_ _type_ _description_ """ + error_ = None load = pd.DataFrame(load_p.sum(axis=1)) total_solar = prod_solar.sum(axis=1) @@ -189,7 +198,15 @@ def generate_economic_dispatch(path_env, start_date_dt, end_date_dt, dt, number_ gens_charac_this["pmax"] = gens_charac_this["Pmax"] gens_charac_this["pmin"] = gens_charac_this["Pmin"] gens_charac_this["cost_per_mw"] = gens_charac_this["marginal_cost"] - economic_dispatch = PypsaDispatcher.from_dataframe(gens_charac_this) + economic_dispatch = PypsaDispatcher.from_dataframe(gens_charac_this, + pypsa_debug_add_inf_gen=pypsa_debug_add_inf_gen, + pypsa_debug_add_thermal=pypsa_debug_add_thermal, + pypsa_debug_add_hydro=pypsa_debug_add_hydro, + pypsa_debug_add_nuc=pypsa_debug_add_nuc, + pypsa_debug_add_solar=pypsa_debug_add_solar, + pypsa_debug_add_wind=pypsa_debug_add_wind, + pypsa_debug_add_slack=pypsa_debug_add_slack, + ) # need to hack it to work... n_gen = len(name_gen) @@ -213,7 +230,7 @@ def generate_economic_dispatch(path_env, start_date_dt, end_date_dt, dt, number_ if res_dispatch is None: error_ = RuntimeError("Pypsa failed to find a solution") - return None, None, None, error_ + return None, None, None, None, error_ # now assign the results final_gen_p = 1.0 * final_gen_p # copy the data frame to avoid modify the original one @@ -232,7 +249,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, hydro_constraints, None + # error_ is None at this stage + return final_gen_p, total_wind_curt, total_solar_curt, hydro_constraints, error_ + def _adjust_gens(all_loss_orig, @@ -255,6 +274,7 @@ def _adjust_gens(all_loss_orig, iter_quality_decrease=50, # acept a reduction of the quality after this number of iteration percentile_quality_decrease=99, hydro_constraints=None, + backend_kwargs_data=None, ): """This function is an auxilliary function. @@ -302,6 +322,8 @@ def _adjust_gens(all_loss_orig, _type_ _description_ """ + if backend_kwargs_data is None: + backend_kwargs_data = {} quality_ = None error_ = None if np.any(~np.isfinite(gen_p)): @@ -354,9 +376,10 @@ def _adjust_gens(all_loss_orig, 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 + assert (target_vector <= 1.0).all(), "possible error detected in the computation of the AC losses: some generators have a target above their pmax!" #### cvxpy - p_t = cp.Variable(shape=(total_step,total_gen)) + p_t = cp.Variable(shape=(total_step, total_gen)) real_p = cp.multiply(p_t, scale_for_loads) constraints = [p_t >= p_min, @@ -409,7 +432,7 @@ def _adjust_gens(all_loss_orig, test=True, # grid_path=grid_path, # assign it the 118 grid param=env_param, - backend=LightSimBackend(), + backend=LightSimBackend(**backend_kwargs_data), chronics_class=FromNPY, # chronics_path=path_chronix2grid, data_feeding_kwargs={"load_p": load_p, @@ -477,17 +500,26 @@ def _adjust_gens(all_loss_orig, def _fix_losses_one_scenario(env_for_loss, - scenario_id, - params, - env_path, - env_param, - load_df, - threshold_stop=0.05, # decide I stop when the data move of less of 0.5 MW at maximum - 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, - ): + scenario_id, + params, + env_path, + env_param, + load_df, + *, + threshold_stop=0.05, # decide I stop when the data move of less of 0.5 MW at maximum + 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, + backend_kwargs_data=None, + pypsa_debug_add_inf_gen=False, + pypsa_debug_add_thermal=True, + pypsa_debug_add_hydro=True, + pypsa_debug_add_nuc=True, + pypsa_debug_add_solar=True, + pypsa_debug_add_wind=True, + pypsa_debug_add_slack=True, + ): """This function is an auxilliary function. Like its main one (see handle_losses) it is here to make sure that if you run an AC model with the data generated, @@ -560,7 +592,15 @@ def _fix_losses_one_scenario(env_for_loss, df["pmax"] = df["Pmax"] df["pmin"] = df["Pmin"] df["cost_per_mw"] = df["marginal_cost"] - economic_dispatch = PypsaDispatcher.from_dataframe(df) + economic_dispatch = PypsaDispatcher.from_dataframe(df, + pypsa_debug_add_inf_gen=pypsa_debug_add_inf_gen, + pypsa_debug_add_thermal=pypsa_debug_add_thermal, + pypsa_debug_add_hydro=pypsa_debug_add_hydro, + pypsa_debug_add_nuc=pypsa_debug_add_nuc, + pypsa_debug_add_solar=pypsa_debug_add_solar, + pypsa_debug_add_wind=pypsa_debug_add_wind, + pypsa_debug_add_slack=pypsa_debug_add_slack, + ) economic_dispatch.read_hydro_guide_curves(os.path.join(ref_pattern_path, 'hydro_french.csv')) economic_dispatch._chronix_scenario = ChroniXScenario(loads=1.0 * load_df, prods=pd.DataFrame(1.0 * gen_p_orig, columns=env_for_loss.name_gen), @@ -598,7 +638,8 @@ def _fix_losses_one_scenario(env_for_loss, max_iter=max_iter, iter_quality_decrease=iter_quality_decrease, percentile_quality_decrease=percentile_quality_decrease, - hydro_constraints=hydro_constraints) + hydro_constraints=hydro_constraints, + backend_kwargs_data=backend_kwargs_data) if error_ is not None: # the procedure failed @@ -607,7 +648,10 @@ 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): + gen_v, start_date_dt, dt_dt, backend_kwargs_data): + if backend_kwargs_data is None: + backend_kwargs_data = {} + with warnings.catch_warnings(): warnings.filterwarnings("ignore") env_for_loss = grid2op.make( @@ -615,21 +659,22 @@ def make_env_for_loss(path_env, env_param, load_p, load_q, final_gen_p, test=True, # grid_path=grid_path, # assign it the 118 grid param=env_param, - backend=LightSimBackend(), + backend=LightSimBackend(**backend_kwargs_data), 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, + "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, + _add_to_name="_env_for_loss" ) return env_for_loss @@ -643,6 +688,7 @@ def handle_losses(path_env, start_date_dt, dt_dt, scenario_id, + *, PmaxErrorCorrRatio=0.9, RampErrorCorrRatio=0.95, threshold_stop=0.05, @@ -650,6 +696,14 @@ def handle_losses(path_env, 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, + backend_kwargs_data=None, + pypsa_debug_add_inf_gen=False, + pypsa_debug_add_thermal=True, + pypsa_debug_add_hydro=True, + pypsa_debug_add_nuc=True, + pypsa_debug_add_solar=True, + pypsa_debug_add_wind=True, + pypsa_debug_add_slack=True, ): """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) @@ -704,11 +758,12 @@ def handle_losses(path_env, 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 ]), + gen_v = np.tile(np.array([float(gens_charac.loc[gens_charac["name"] == nm_gen].V.iloc[0]) 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) + final_gen_p, gen_v, start_date_dt, dt_dt, + backend_kwargs_data=backend_kwargs_data) res_gen_p, error_, quality_ = _fix_losses_one_scenario(env_for_loss, scenario_id, loss_param, @@ -722,6 +777,14 @@ def handle_losses(path_env, # replace the "at maximum" by "percentile 99%" percentile_quality_decrease=percentile_quality_decrease, hydro_constraints=hydro_constraints, + backend_kwargs_data=backend_kwargs_data, + pypsa_debug_add_inf_gen=pypsa_debug_add_inf_gen, + pypsa_debug_add_thermal=pypsa_debug_add_thermal, + pypsa_debug_add_hydro=pypsa_debug_add_hydro, + pypsa_debug_add_nuc=pypsa_debug_add_nuc, + pypsa_debug_add_solar=pypsa_debug_add_solar, + pypsa_debug_add_wind=pypsa_debug_add_wind, + pypsa_debug_add_slack=pypsa_debug_add_slack, ) if error_ is not None: return None, error_, None, env_for_loss @@ -902,6 +965,7 @@ def save_meta_data(this_scen_path, 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, gen_type, @@ -912,6 +976,7 @@ def generate_a_scenario(path_env, load_seed, renew_seed, gen_p_forecast_seed, + *, handle_loss=True, nb_steps=None, PmaxErrorCorrRatio=0.9, @@ -924,6 +989,14 @@ def generate_a_scenario(path_env, tol_zero=1e-3, debug=True, # TODO more feature ! load_weekly_pattern=None, + backend_kwargs_data=None, + pypsa_debug_add_inf_gen=False, + pypsa_debug_add_thermal=True, + pypsa_debug_add_hydro=True, + pypsa_debug_add_nuc=True, + pypsa_debug_add_solar=True, + pypsa_debug_add_wind=True, + pypsa_debug_add_slack=True, ): """This function generates and save the data for a scenario. @@ -1014,7 +1087,6 @@ def generate_a_scenario(path_env, 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 @@ -1043,13 +1115,32 @@ def generate_a_scenario(path_env, 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) + final_gen_p, gens_charac, opf_params, + pypsa_debug_add_inf_gen=pypsa_debug_add_inf_gen, + pypsa_debug_add_thermal=pypsa_debug_add_thermal, + pypsa_debug_add_hydro=pypsa_debug_add_hydro, + pypsa_debug_add_nuc=pypsa_debug_add_nuc, + pypsa_debug_add_solar=pypsa_debug_add_solar, + pypsa_debug_add_wind=pypsa_debug_add_wind, + pypsa_debug_add_slack=pypsa_debug_add_slack, + ) 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 + if ((pypsa_debug_add_inf_gen is not False) or + (pypsa_debug_add_thermal is not True) or + (pypsa_debug_add_hydro is not True) or + (pypsa_debug_add_nuc is not True) or + (pypsa_debug_add_solar is not True) or + (pypsa_debug_add_wind is not True) or + (pypsa_debug_add_slack is not True)): + # one of the pypsa flag is not to its default, data cannot be generated like that + error_ = ("One of the pypsa flag is set to a 'debug' mode, data will not be generated further. " + "If you have a problem, it should not be related to the 'OPF'") + 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: @@ -1063,11 +1154,19 @@ def generate_a_scenario(path_env, start_date_dt, dt_dt, scenario_id, + pypsa_debug_add_inf_gen=pypsa_debug_add_inf_gen, + pypsa_debug_add_thermal=pypsa_debug_add_thermal, + pypsa_debug_add_hydro=pypsa_debug_add_hydro, + pypsa_debug_add_nuc=pypsa_debug_add_nuc, + pypsa_debug_add_solar=pypsa_debug_add_solar, + pypsa_debug_add_wind=pypsa_debug_add_wind, + pypsa_debug_add_slack=pypsa_debug_add_slack, PmaxErrorCorrRatio=PmaxErrorCorrRatio, RampErrorCorrRatio=RampErrorCorrRatio, threshold_stop=threshold_stop, max_iter=max_iter, - hydro_constraints=hydro_constraints) + hydro_constraints=hydro_constraints, + backend_kwargs_data=backend_kwargs_data) if error_ is not None: # TODO log that ! return error_, None, None, None, None, None, None, None @@ -1082,7 +1181,8 @@ def generate_a_scenario(path_env, 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) + start_date_dt, dt_dt, + backend_kwargs_data=backend_kwargs_data) prng = default_rng(gen_p_forecast_seed) beg_forca = time.perf_counter() From fc5bc341cdacedb80ee02773174afe38e6a7ebba Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Fri, 5 Apr 2024 18:08:59 +0200 Subject: [PATCH 3/8] fixing some issue and improve 'convergence' of losses --- .../run_economic_dispatch.py | 14 +- .../_EDispatch_L2RPN2020/utils.py | 23 ++- chronix2grid/grid2op_utils/utils.py | 162 +++++++++++------- 3 files changed, 126 insertions(+), 73 deletions(-) diff --git a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/run_economic_dispatch.py b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/run_economic_dispatch.py index 455dc4da..ba6f6280 100644 --- a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/run_economic_dispatch.py +++ b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/run_economic_dispatch.py @@ -11,6 +11,7 @@ import time from typing import Dict, Literal +import numpy as np import pandas as pd import pypsa @@ -95,8 +96,10 @@ def main_run_disptach(pypsa_net, # user specified a pmin for the slack bus slack_pmin = float(params["slack_pmin"]) / float(pypsa_net.generators.loc[slack_name].p_nom) if slack_name not in gen_constraints["p_min_pu"]: - print(gen_constraints["p_min_pu"].columns) - raise RuntimeError("Unknown error in chronix2grid: the provided slack has no constraints on the pypsa grid...") + gen_constraints["p_min_pu"][slack_name] = slack_pmin + else: + gen_constraints["p_min_pu"][slack_name] = np.maximum(slack_pmin, gen_constraints["p_min_pu"][slack_name].values) + if slack_pmin < gen_constraints["p_min_pu"][slack_name].min(): raise RuntimeError("Wrong configuration detected: you put in 'params_opf' " "a minimum value for the slack ('slack_pmin') below its pmin.") @@ -106,14 +109,15 @@ def main_run_disptach(pypsa_net, # user specified a pmin for the slack bus slack_pmax = float(params["slack_pmax"]) / float(pypsa_net.generators.loc[slack_name].p_nom) if slack_name not in gen_constraints["p_max_pu"]: - print(gen_constraints["p_min_pu"].columns) - raise RuntimeError("Unknown error in chronix2grid: the provided slack has no constraints on the pypsa grid...") + gen_constraints["p_max_pu"][slack_name] = slack_pmax + else: + gen_constraints["p_max_pu"][slack_name] = np.minimum(slack_pmin, gen_constraints["p_max_pu"][slack_name].values) if slack_pmax > gen_constraints["p_max_pu"][slack_name].max(): raise RuntimeError("Wrong configuration detected: you put in 'params_opf' " "a maximum value for the slack ('slack_pmax') above its pmax.") if slack_pmax > 1.: print("Doubtful parameters: slack pmax would be higher than its pmax, basically " - "the parameters `slack_pmin` of `params_opf.json` is ignored.") + "the parameters `slack_pmax` of `params_opf.json` is ignored.") error_ = False start = time.time() diff --git a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py index cfb209aa..82ab88cb 100644 --- a/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py +++ b/chronix2grid/generation/_dispatch/_PypsaDispatchBackend/_EDispatch_L2RPN2020/utils.py @@ -424,19 +424,26 @@ def run_opf(net, # raise error for non feasible configurations if (net.generators["p_min_pu"].values < 0.).any(): - raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative p_min for {net.generators.loc[net.generators['p_min_pu'].values < 0.].index}") + raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative " + f"p_min for {net.generators.loc[net.generators['p_min_pu'].values < 0.].index}") if (net.generators["p_max_pu"].values < 0.).any(): - raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative p_max for {net.generators.loc[net.generators['p_max_pu'].values < 0.].index}") + raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative " + f"p_max for {net.generators.loc[net.generators['p_max_pu'].values < 0.].index}") if (net.generators["ramp_limit_up"].values < 0.).any(): - raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative ramps (up) for {net.generators.loc[net.generators['ramp_limit_up'].values < 0.].index}") + raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative " + f"ramps (up) for {net.generators.loc[net.generators['ramp_limit_up'].values < 0.].index}") if (net.generators["ramp_limit_down"].values < 0.).any(): - raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative ramps (down) for {net.generators.loc[net.generators['ramp_limit_down'].values < 0.].index}") + raise RuntimeError(f"Wrong input data (prods_charac.csv) there are some negative ramps " + f"(down) for {net.generators.loc[net.generators['ramp_limit_down'].values < 0.].index}") if (net.generators_t.p_min_pu.min() > 1.).any(): debug_info = net.generators_t.p_min_pu.min() - raise RuntimeError(f"A generator would have a p_min above its pmax ! most likely the slack... Check {debug_info.loc[debug_info > 1.].index.values}") + raise RuntimeError(f"A generator would have a p_min above its pmax ! most likely " + f"the slack... Check {debug_info.loc[debug_info > 1.].index.values}") if ((net.generators_t["p_max_pu"] - net.generators_t["p_min_pu"]) < 0.).any().any(): tmp_debug = ((net.generators_t["p_max_pu"] - net.generators_t["p_min_pu"]) < 0.).any() - raise RuntimeError(f"Invalid configuration: infeasibility detected for some generators (most likely the slack), check generator {tmp_debug[tmp_debug].index.values}") + raise RuntimeError(f"Invalid configuration: infeasibility detected for some " + f"generators (most likely the slack), check generator " + f"{tmp_debug[tmp_debug].index.values}") # try to detect some non feasible states if "agg_wind" in net.generators_t["p_max_pu"]: @@ -454,7 +461,9 @@ def run_opf(net, if (min_net_demand > total_pmax).any(): return None, RuntimeError("Some non feasible time steps are found (net demand above pmax of controlable generators)") if (demand["agg_load"].values < total_pmin).any(): - return None, RuntimeError("Some non feasible time steps are found (demand [with full curtailment of solar and wind] is below pmin of controlable generators)") + return None, RuntimeError("Some non feasible time steps are found (demand [with full " + "curtailment of solar and wind] is below pmin of " + "controlable generators)") # TODO check also feasible ramps # tmp_ = net.generators.loc[(net.generators.index != "agg_wind") & (net.generators.index != "agg_solar"), ["p_nom", "ramp_limit_up", "ramp_limit_down"]].values diff --git a/chronix2grid/grid2op_utils/utils.py b/chronix2grid/grid2op_utils/utils.py index 07b7885e..28c461c6 100644 --- a/chronix2grid/grid2op_utils/utils.py +++ b/chronix2grid/grid2op_utils/utils.py @@ -13,6 +13,7 @@ import json import shutil import time +import grid2op.Environment import pandas as pd import os import cvxpy as cp @@ -23,6 +24,7 @@ from grid2op.Chronics import FromNPY from grid2op.Action import DontAct from grid2op.Opponent import NeverAttackBudget, BaseOpponent +from grid2op.Exceptions import Grid2OpException from lightsim2grid import LightSimBackend @@ -388,15 +390,28 @@ def _adjust_gens(all_loss_orig, p_t[1:,:] - p_t[:-1,:] <= ramp_max, cp.sum(real_p, axis=1) == load.reshape(-1), ] - cost = cp.sum_squares(p_t - target_vector) + cp.norm1(cp.multiply(p_t, turned_off_orig)) + # for some environment (on more realistic data), we need to put a constant before the norm1 + # otherwise it diverges + cost = cp.sum_squares(p_t - target_vector) + 1e-2 * cp.norm1(cp.multiply(p_t, turned_off_orig)) prob = cp.Problem(cp.Minimize(cost), constraints) try: 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 + # we try different tolerance to increase the "likelyhood" + # of convergence, and also because we do not really care in this routine + # to be optimal + print("cvxpy failed with given tolerance, trying to relax them to 1e-4") + try: + res = prob.solve(verbose=True, eps_abs=1e-4, eps_rel=1e-4) + except cp.error.SolverError as exc_: + print("cvxpy failed with given tolerance, trying to relax them to 1e-3") + try: + res = prob.solve(verbose=True, eps_abs=1e-3, eps_rel=1e-3) + 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") @@ -425,26 +440,33 @@ def _adjust_gens(all_loss_orig, id_redisp += 1 # 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(**backend_kwargs_data), - 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}, - opponent_budget_per_ts=0., - opponent_init_budget=0., - opponent_class=BaseOpponent, - opponent_budget_class=NeverAttackBudget, - opponent_action_class=DontAct, - ) + try: + 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(**backend_kwargs_data), + 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}, + opponent_budget_per_ts=0., + opponent_init_budget=0., + opponent_class=BaseOpponent, + opponent_budget_class=NeverAttackBudget, + opponent_action_class=DontAct, + ) + except Grid2OpException as exc_: + error_ = RuntimeError(f"grid2op cannot create the environment at iteration {iter_num}, error {exc_}") + res_gen_p = None + quality_ = None + break + diff_ = np.full((env_fixed.max_episode_duration(), env_fixed.n_gen), fill_value=np.NaN) all_loss[:] = np.NaN @@ -458,6 +480,12 @@ def _adjust_gens(all_loss_orig, obs, reward, done, info = env_fixed.step(env_fixed.action_space()) i += 1 if done: + if info["exception"]: + # there is a powerflow error here + res_gen_p = None + quality_ = None + error_ = RuntimeError(f"lightsim2grid cannot compute powerflow at iteration {iter_num}, error {info['exception']}") + return res_gen_p, error_, quality_ break all_loss[i] = np.sum(obs.gen_p) - np.sum(obs.load_p) diff_[i] = obs.gen_p - res_gen_p[i] @@ -499,7 +527,7 @@ def _adjust_gens(all_loss_orig, return res_gen_p, error_, quality_ -def _fix_losses_one_scenario(env_for_loss, +def _fix_losses_one_scenario(env_for_loss : grid2op.Environment.Environment, scenario_id, params, env_path, @@ -574,6 +602,9 @@ def _fix_losses_one_scenario(env_for_loss, while not done: obs, reward, done, info = env_for_loss.step(env_for_loss.action_space()) if done: + if info["exception"]: + # there is a powerflow error here + return None, f"Error at first losses computation: {info['exception']}", None break i += 1 all_loss_orig[i] = np.sum(obs.gen_p) - np.sum(obs.load_p) @@ -651,32 +682,35 @@ def make_env_for_loss(path_env, env_param, load_p, load_q, final_gen_p, gen_v, start_date_dt, dt_dt, backend_kwargs_data): if backend_kwargs_data is None: backend_kwargs_data = {} - - 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(**backend_kwargs_data), - 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, - _add_to_name="_env_for_loss" - ) - return env_for_loss + try: + 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(**backend_kwargs_data), + 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, + _add_to_name="_env_for_loss" + ) + except Grid2OpException as exc_: + return None, exc_ + return env_for_loss, None + def handle_losses(path_env, n_gen, @@ -761,9 +795,13 @@ def handle_losses(path_env, gen_v = np.tile(np.array([float(gens_charac.loc[gens_charac["name"] == nm_gen].V.iloc[0]) 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, - backend_kwargs_data=backend_kwargs_data) + env_for_loss, exc_ = make_env_for_loss(path_env, env_param, load_p, load_q, + final_gen_p, gen_v, start_date_dt, dt_dt, + backend_kwargs_data=backend_kwargs_data) + + if exc_ is not None: + return None, exc_, None, None + res_gen_p, error_, quality_ = _fix_losses_one_scenario(env_for_loss, scenario_id, loss_param, @@ -1178,12 +1216,14 @@ def generate_a_scenario(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) - 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, - backend_kwargs_data=backend_kwargs_data) - + env_for_loss, exc_ = make_env_for_loss(path_env, env_param, + load_p, load_q, + final_gen_p, gen_v, + start_date_dt, dt_dt, + backend_kwargs_data=backend_kwargs_data) + if exc_ is not None: + return exc_, None, None, None, None, None, None, None + prng = default_rng(gen_p_forecast_seed) beg_forca = time.perf_counter() tmp_ = generate_forecasts_gen(new_forecasts, From e0b08df1a1a915a5e8751578281a5c458baa56c8 Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Tue, 30 Jul 2024 11:24:21 +0200 Subject: [PATCH 4/8] adding tests for grid2op integration --- .circleci/config.yml | 4 + .../grid2op_utils/generate_one_episode.py | 14 ++- tests/grid2op_tests/__init__.py | 0 .../grid2op_tests/test_grid2op_integration.py | 119 ++++++++++++++++++ 4 files changed, 135 insertions(+), 2 deletions(-) create mode 100644 tests/grid2op_tests/__init__.py create mode 100644 tests/grid2op_tests/test_grid2op_integration.py diff --git a/.circleci/config.yml b/.circleci/config.yml index 4439f7cd..33f69fd6 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -46,6 +46,10 @@ jobs: name: Run integration tests # This assumes pytest is installed via the install-package step above command: pytest tests/integration_tests + - run: + name: Run some grid2op tests (test that it does not break grid2op `FromChronix2grid`) + # This assumes pytest is installed via the install-package step above + command: pytest tests/grid2op_tests # Invoke jobs via workflows diff --git a/chronix2grid/grid2op_utils/generate_one_episode.py b/chronix2grid/grid2op_utils/generate_one_episode.py index 0d6e1f77..fb33e09b 100644 --- a/chronix2grid/grid2op_utils/generate_one_episode.py +++ b/chronix2grid/grid2op_utils/generate_one_episode.py @@ -57,8 +57,18 @@ def generate_one_episode(env: grid2op.Environment.Environment, error_ = "first" 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, + res_gen = generate_a_scenario(path_env=path_env, + name_gen=name_gen, + gen_type=gen_type, + output_dir=output_dir, + start_date=start_date, + dt=dt, + scen_id=scen_id, + load_seed=load_seed, + renew_seed=renew_seed, + gen_p_forecast_seed=gen_p_forecast_seed, + handle_loss=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/tests/grid2op_tests/__init__.py b/tests/grid2op_tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/grid2op_tests/test_grid2op_integration.py b/tests/grid2op_tests/test_grid2op_integration.py new file mode 100644 index 00000000..5b6a9422 --- /dev/null +++ b/tests/grid2op_tests/test_grid2op_integration.py @@ -0,0 +1,119 @@ +# 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 Grid2Op, Grid2Op a testbed platform to model sequential decision making in power systems. + +import pdb +import warnings +import os +import grid2op +import numpy as np +from grid2op.Chronics import FromChronix2grid +import unittest +import pkg_resources +from lightsim2grid import LightSimBackend + +DEV_DATA_FOLDER = pkg_resources.resource_filename("grid2op", "data") + +class TestFromChronix2Grid(unittest.TestCase): + def _aux_reset_env(self): + self.env.seed(self.seed_) + self.env.set_id(self.scen_nm) + return self.env.reset() + + def setUp(self) -> None: + self.seed_ = 0 + self.env_nm = "l2rpn_wcci_2022_dev" + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + self.env = grid2op.make(self.env_nm, + test=True, + backend=LightSimBackend(), + chronics_class=FromChronix2grid, + data_feeding_kwargs={"env_path": os.path.join(DEV_DATA_FOLDER, self.env_nm), # os.path.join(grid2op.get_current_local_dir(), self.env_nm), + "with_maintenance": True, + "max_iter": 10, + "with_loss": False + }, + _add_to_name=type(self).__name__ + ) + + + def test_ok(self): + """test it can be created""" + assert self.env.chronics_handler.real_data._gen_p.shape == (12, 62) + assert np.all(np.isfinite(self.env.chronics_handler.real_data._gen_p)) + + def test_seed_setid(self): + """test env.seed(...) and env.set_id(...)""" + id_ref = '2525122259@2050-02-28' + id_ref = '377638611@2050-02-28' + # test tell_id + sum_prod_ref = 42340.949878 + sum_prod_ref = 41784.477161 + self.env.seed(self.seed_) + self.env.reset() + id_ = self.env.chronics_handler.get_id() + assert id_ == id_ref, f"wrong id {id_} instead of {id_ref}" + assert abs(self.env.chronics_handler.real_data._gen_p.sum() - sum_prod_ref) <= 1e-4, f"{self.env.chronics_handler.real_data._gen_p.sum():.2f}" + self.env.reset() + # assert abs(self.env.chronics_handler.real_data._gen_p.sum() - 38160.833356999996) <= 1e-4 + assert abs(self.env.chronics_handler.real_data._gen_p.sum() - 37662.206248999995) <= 1e-4, f"{self.env.chronics_handler.real_data._gen_p.sum():.2f}" + self.env.set_id(id_ref) + self.env.reset() + assert abs(self.env.chronics_handler.real_data._gen_p.sum() - sum_prod_ref) <= 1e-4, f"{self.env.chronics_handler.real_data._gen_p.sum():.2f}" + + # test seed + self.env.seed(self.seed_) + self.env.reset() + assert abs(self.env.chronics_handler.real_data._gen_p.sum() - sum_prod_ref) <= 1e-4, f"{self.env.chronics_handler.real_data._gen_p.sum():.2f}" + + def test_episode(self): + """test that the episode can go until the end""" + self.env.seed(1) + obs = self.env.reset() + assert obs.current_step == 0 + assert obs.max_step == 10 + for i in range(obs.max_step - 1): + obs, reward, done, info = self.env.step(self.env.action_space()) + assert not done, f"episode should not be over at iteration {i}" + obs, reward, done, info = self.env.step(self.env.action_space()) + assert obs.current_step == 10 + assert obs.max_step == 10 + assert done, "episode should have terminated" + obs = self.env.reset() + assert obs.max_step == 10 + assert obs.current_step == 0 + + def test_maintenance(self): + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + self.env = grid2op.make(self.env_nm, + test=True, + _add_to_name=type(self).__name__, + chronics_class=FromChronix2grid, + data_feeding_kwargs={"env_path": os.path.join(DEV_DATA_FOLDER, self.env_nm), + "with_maintenance": True, + "max_iter": 2 * 288, + "with_loss": False + } + ) + self.env.seed(0) + id_ref = '0@2050-08-08' + self.env.set_id(id_ref) + obs = self.env.reset() + assert np.all(obs.time_next_maintenance[[43, 126]] == [107, 395]) + assert np.all(obs.time_next_maintenance[:43] == -1) + assert np.all(obs.time_next_maintenance[127:] == -1) + assert np.all(obs.time_next_maintenance[44:126] == -1) + assert self.env.chronics_handler.real_data.maintenance is not None + assert self.env.chronics_handler.real_data.maintenance.sum() == 192 + assert self.env.chronics_handler.real_data.maintenance_time is not None + assert self.env.chronics_handler.real_data.maintenance_duration is not None + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file From f455ae1145b1bc61b318c6b072a6781e0f808855 Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Tue, 30 Jul 2024 16:20:39 +0200 Subject: [PATCH 5/8] start to add feature to have one pattern per loads [skip ci] --- .../consumption/consumption_utils.py | 153 ++++++++++++++---- .../generation/consumption/generate_load.py | 9 +- chronix2grid/grid2op_utils/loads_utils.py | 20 ++- 3 files changed, 144 insertions(+), 38 deletions(-) diff --git a/chronix2grid/generation/consumption/consumption_utils.py b/chronix2grid/generation/consumption/consumption_utils.py index 2acbc504..24602a59 100755 --- a/chronix2grid/generation/consumption/consumption_utils.py +++ b/chronix2grid/generation/consumption/consumption_utils.py @@ -14,19 +14,40 @@ 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, +__CACHED__DATA = {} # DO NOT MODIFY THIS UNDER ANY CIRCUMSTANCES ! +# TODO optim: cache the data instead of re reading everything every time + + +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 isinstance(load_weekly_pattern, pd.DataFrame): + weekly_pattern_vals = load_weekly_pattern['test'].values + weekly_pattern_datetimes = load_weekly_pattern["datetime"] + elif isinstance(load_weekly_pattern, os.PathLike): + weekly_pattern_vals = load_weekly_pattern + if not (load_weekly_pattern / "datetimes.csv").exists(): + raise RuntimeError(f"Impossible to load the associated dates and times in {load_weekly_pattern}. " + f"Make sure that the files `datetimes.csv` exists there (and is a pandas " + f"dataframe with at least a column name 'datetime')") + weekly_pattern_datetimes = pd.read_csv(load_weekly_pattern / "datetimes.csv", sep=";")["datetime"] + else: + raise RuntimeError("at this stage `load_weekly_pattern` args should be either a path or a dataframe") + 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") + datetime_lwp = pd.to_datetime(weekly_pattern_datetimes, 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]) @@ -37,23 +58,32 @@ def compute_loads(loads_charac, temperature_noise, params, load_weekly_pattern, # day_lag = 0 loads_series = {} ref_curves = None - for i, name in enumerate(loads_charac['name']): - mask = (loads_charac['name'] == name) + for i, load_name in enumerate(loads_charac['name']): + mask = (loads_charac['name'] == load_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] - tmp_ = 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_vals, + index=i, + day_lag=day_lag, + add_dim=add_dim, return_ref_curve=return_ref_curve, isoweekday_lwp=isoweekday, - hour_minutes_lwp=hour_minutes) + hour_minutes_lwp=hour_minutes, + load_name=load_name, + weekly_pattern_datetimes=weekly_pattern_datetimes) 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_ + loads_series[load_name], ref_curves[:,i] = tmp_ else: - loads_series[name] = tmp_ + loads_series[load_name] = tmp_ + else: + raise NotImplementedError(f"Load generation is not implemented for type {loads_charac[mask]['type'].values}") if loads_charac[mask]['type'].values == 'industrial': raise NotImplementedError("Impossible to generate industrial loads for now.") @@ -77,11 +107,19 @@ def get_seasonal_pattern(params): 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): + hour_minutes_lwp=None, + load_name=None, + weekly_pattern_datetimes=None): # Compute refined signals @@ -97,7 +135,13 @@ def compute_residential(locations, Pmax, temperature_noise, params, seasonal_pattern = get_seasonal_pattern(params) # Get weekly pattern - weekly_pattern = compute_load_pattern(params, weekly_pattern, index, day_lag, isoweekday_lwp=isoweekday_lwp, hour_minutes_lwp=hour_minutes_lwp) + weekly_pattern = compute_load_pattern(params, + weekly_pattern, + index, day_lag, + isoweekday_lwp=isoweekday_lwp, + hour_minutes_lwp=hour_minutes_lwp, + load_name=load_name, + weekly_pattern_datetimes=weekly_pattern_datetimes) std_temperature_noise = params['std_temperature_noise'] residential_series = Pmax * weekly_pattern * (std_temperature_noise * temperature_signal + seasonal_pattern) @@ -105,7 +149,14 @@ def compute_residential(locations, Pmax, temperature_noise, params, return residential_series, Pmax * weekly_pattern * seasonal_pattern return residential_series -def compute_load_pattern(params, weekly_pattern, index, day_lag=None, isoweekday_lwp=None, hour_minutes_lwp=None): +def compute_load_pattern(params, + weekly_pattern, + index, + day_lag=None, + isoweekday_lwp=None, + hour_minutes_lwp=None, + load_name=None, + weekly_pattern_datetimes=None): """ Loads a typical hourly pattern, and interpolates it to generate a smooth solar generation pattern between 0 and 1 @@ -118,30 +169,69 @@ def compute_load_pattern(params, weekly_pattern, index, day_lag=None, isoweekday Output: (np.array) A smooth solar pattern """ + index_weekly_perweek = 12 * 24 * 7 + generate_multiple_pattern = None # solar_pattern resolution : 1H, 8761 - + if isinstance(weekly_pattern, os.PathLike): + # user provided a pattern in the shape of a folder + if load_name is None: + raise RuntimeError("When you want to use a different load pattern for each load, " + "you should provide the load name, ie the args `load_name`.") + if weekly_pattern_datetimes is None: + raise RuntimeError("When you want to use a different load pattern for each load, " + "you should provide the dates and times for this reference data " + "ie the args `weekly_pattern_datetimes`.") + + pattern_path = weekly_pattern / f"{load_name}.csv" + if not pattern_path.exists(): + pattern_path = weekly_pattern / f"{load_name}.csv.gz" + if not pattern_path.exists(): + raise RuntimeError(f"Unable to locate the weekly pattern for load {load_name} at {weekly_pattern}") + tmp_weekly_pattern = pd.read_csv(pattern_path, sep=";") + used_weekly_pattern = 1. * tmp_weekly_pattern["values"].values + if isoweekday_lwp is None or hour_minutes_lwp is None: + raise RuntimeError(f"When provided a different pattern for each loads, you should also provide " + f"the isoweekday to simulate and also the hours / minute") + generate_multiple_pattern = False + elif isinstance(weekly_pattern, np.ndarray): + # standard legacy behaviour + used_weekly_pattern = 1. * weekly_pattern + generate_multiple_pattern = True + else: + raise RuntimeError(f"weekly_pattern should be either a path or a numpy array, you provided {type(weekly_pattern)}") + # Keep only one week of pattern - index_weekly_perweek = 12 * 24 * 7 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) + index %= int((used_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] - + isoweekday_start = params["start_date"].isoweekday() + if generate_multiple_pattern: + # in this case only one pattern is given for all the loads, so I come up with a solution to generate + # different patterns from this one + 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] + else: + # one pattern per load is already provided so I only here to extract the correct "week" + week_of_year_params = params["start_date"].isocalendar()[1] # the week number of the generated data + week_of_year_data = np.array([el.week for el in pd.to_datetime(weekly_pattern_datetimes, format="%Y-%m-%d %H:%M:%S")]) + possible_first_index = np.where((week_of_year_data == week_of_year_params) & + (isoweekday_lwp == isoweekday_start) & + (iso_hm_start == hour_minutes_lwp))[0] + first_index = possible_first_index[0] + # 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) + used_weekly_pattern = 1.0 * used_weekly_pattern[first_index:last_index] + used_weekly_pattern /= np.mean(used_weekly_pattern) # now generate an ouput of the right length if isoweekday_lwp is None or hour_minutes_lwp is None: @@ -150,8 +240,8 @@ def compute_load_pattern(params, weekly_pattern, index, day_lag=None, isoweekday 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 = np.tile(weekly_pattern, N_repet) + N_repet = int((Nt_inter_hr - 1) // len(used_weekly_pattern) + 1) + stacked_weekly_pattern = np.tile(used_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) @@ -167,10 +257,9 @@ def compute_load_pattern(params, weekly_pattern, index, day_lag=None, isoweekday 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) + N_repet = np.ceil(nb_ts / used_weekly_pattern.shape[0]).astype(int) + stacked_weekly_pattern = np.tile(used_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 e7ce0b13..e8c0aad3 100755 --- a/chronix2grid/generation/consumption/generate_load.py +++ b/chronix2grid/generation/consumption/generate_load.py @@ -31,8 +31,13 @@ def get_add_dim(params, loads_charac): return add_dim -def main(scenario_destination_path, seed, params, loads_charac, - load_weekly_pattern, write_results = True, day_lag=0, +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): """ diff --git a/chronix2grid/grid2op_utils/loads_utils.py b/chronix2grid/grid2op_utils/loads_utils.py index 3c0dc32e..748e544e 100644 --- a/chronix2grid/grid2op_utils/loads_utils.py +++ b/chronix2grid/grid2op_utils/loads_utils.py @@ -10,11 +10,12 @@ import json import pandas as pd import numpy as np +import pathlib 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) + compute_load_pattern) from chronix2grid.generation.consumption.generate_load import get_add_dim as get_add_dim_load @@ -71,7 +72,10 @@ def generate_new_loads(load_seed, 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") + if isinstance(load_weekly_pattern, pd.DataFrame): + datetime_lwp = pd.to_datetime(load_weekly_pattern["datetime"], format="%Y-%m-%d %H:%M:%S") + else: + raise NotImplementedError("TODO load weekly pattern !") 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, @@ -181,9 +185,17 @@ def generate_loads(path_env, if load_weekly_pattern is None: load_weekly_pattern = pd.read_csv(os.path.join(ref_pattern_path, "load_weekly_pattern.csv"), sep=",") else: - load_weekly_pattern = pd.DataFrame(load_weekly_pattern) + if isinstance(load_weekly_pattern, (str, os.PathLike)): + load_weekly_pattern = pathlib.Path(load_weekly_pattern) + elif isinstance(load_weekly_pattern, (np.ndarray, pd.DataFrame)): + load_weekly_pattern = pd.DataFrame(load_weekly_pattern) + else: + raise TypeError("load_weekly_pattern kwargs should be either a str " + "(or pathlike) or a pandas Dataframe (or numpy array) and not " + f"{type(load_weekly_pattern)}") if new_forecasts: + # TODO load weekly pattern ! load_p, load_p_forecasted, load_ref_curve = generate_new_loads(load_seed, load_params, forecasts_params, @@ -201,7 +213,7 @@ def generate_loads(path_env, 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, + load_p, load_p_forecasted, load_ref_curve = load_generator.run(load_weekly_pattern=load_weekly_pattern, # TODO load weekly pattern return_ref_curve=True, use_legacy=False) From 0ce744496fdb73f7c5334c23713a05c1a03c3122 Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Fri, 2 Aug 2024 16:53:23 +0200 Subject: [PATCH 6/8] chronix2grid able to digest individual load pattern --- .../consumption/consumption_utils.py | 38 +++++++++++++++++-- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/chronix2grid/generation/consumption/consumption_utils.py b/chronix2grid/generation/consumption/consumption_utils.py index 24602a59..171f37dd 100755 --- a/chronix2grid/generation/consumption/consumption_utils.py +++ b/chronix2grid/generation/consumption/consumption_utils.py @@ -156,7 +156,10 @@ def compute_load_pattern(params, isoweekday_lwp=None, hour_minutes_lwp=None, load_name=None, - weekly_pattern_datetimes=None): + weekly_pattern_datetimes=None, + use_cache=True, # TODO propagate to function calling this one + fill_na=True, # TODO propagate to function calling this one + ): """ Loads a typical hourly pattern, and interpolates it to generate a smooth solar generation pattern between 0 and 1 @@ -187,8 +190,21 @@ def compute_load_pattern(params, pattern_path = weekly_pattern / f"{load_name}.csv.gz" if not pattern_path.exists(): raise RuntimeError(f"Unable to locate the weekly pattern for load {load_name} at {weekly_pattern}") - tmp_weekly_pattern = pd.read_csv(pattern_path, sep=";") - used_weekly_pattern = 1. * tmp_weekly_pattern["values"].values + if not use_cache: + tmp_weekly_pattern = pd.read_csv(pattern_path, sep=";") + if tmp_weekly_pattern.isna().any().any(): + print(f"WARNING: NA detected for load pattern {pattern_path}") + else: + if pattern_path in __CACHED__DATA: + tmp_weekly_pattern = __CACHED__DATA[pattern_path] + print(f"Using cached data for {pattern_path}") + else: + tmp_weekly_pattern = pd.read_csv(pattern_path, sep=";") + if tmp_weekly_pattern.isna().any().any(): + print(f"WARNING: NA detected for load pattern {pattern_path}") + print(f"Loading cached data for {pattern_path}") + __CACHED__DATA[pattern_path] = tmp_weekly_pattern + used_weekly_pattern = tmp_weekly_pattern["values"].values if isoweekday_lwp is None or hour_minutes_lwp is None: raise RuntimeError(f"When provided a different pattern for each loads, you should also provide " f"the isoweekday to simulate and also the hours / minute") @@ -231,7 +247,21 @@ def compute_load_pattern(params, # now extract right week of data last_index = first_index + index_weekly_perweek used_weekly_pattern = 1.0 * used_weekly_pattern[first_index:last_index] - used_weekly_pattern /= np.mean(used_weekly_pattern) + not_finite = ~np.isfinite(used_weekly_pattern) + if not_finite.any(): + print(f"WARNING: NA detected for load pattern for load index {index} (name {load_name})") + if fill_na: + if not not_finite.all(): + used_weekly_pattern[not_finite] = np.nanmean(used_weekly_pattern) + else: + print(f"WARNING: pattern full of nan for for load index {index} (name {load_name}), replaced by 1.") + used_weekly_pattern[:] = 1. + avg_ = np.mean(used_weekly_pattern) + if np.abs(avg_) > 1e-5: + used_weekly_pattern /= avg_ + else: + print(f"WARNING: for load index {index} (name {load_name}): average value too small, pattern replaced by all one") + used_weekly_pattern[:] = 1. # now generate an ouput of the right length if isoweekday_lwp is None or hour_minutes_lwp is None: From 07104987b09aaa0a758dd1a9046dc2c59a46e558 Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Fri, 8 Nov 2024 10:53:20 +0100 Subject: [PATCH 7/8] try to adress #83 --- chronix2grid/generation/consumption/consumption_utils.py | 2 -- chronix2grid/generation/renewable/solar_wind_utils.py | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/chronix2grid/generation/consumption/consumption_utils.py b/chronix2grid/generation/consumption/consumption_utils.py index 171f37dd..bd73177f 100755 --- a/chronix2grid/generation/consumption/consumption_utils.py +++ b/chronix2grid/generation/consumption/consumption_utils.py @@ -197,12 +197,10 @@ def compute_load_pattern(params, else: if pattern_path in __CACHED__DATA: tmp_weekly_pattern = __CACHED__DATA[pattern_path] - print(f"Using cached data for {pattern_path}") else: tmp_weekly_pattern = pd.read_csv(pattern_path, sep=";") if tmp_weekly_pattern.isna().any().any(): print(f"WARNING: NA detected for load pattern {pattern_path}") - print(f"Loading cached data for {pattern_path}") __CACHED__DATA[pattern_path] = tmp_weekly_pattern used_weekly_pattern = tmp_weekly_pattern["values"].values if isoweekday_lwp is None or hour_minutes_lwp is None: diff --git a/chronix2grid/generation/renewable/solar_wind_utils.py b/chronix2grid/generation/renewable/solar_wind_utils.py index 8ea7a137..443a4f60 100755 --- a/chronix2grid/generation/renewable/solar_wind_utils.py +++ b/chronix2grid/generation/renewable/solar_wind_utils.py @@ -126,7 +126,7 @@ def compute_solar_pattern(params, solar_pattern, tol=0.0): end_min = int(pd.Timedelta(params['end_date'] - start_year).total_seconds() // 60) Nt_inter_hr = int(end_min // 60 + 1) - N_repet = int((Nt_inter_hr - 1) // len(solar_pattern) + 1) + N_repet = int((Nt_inter_hr) // len(solar_pattern) + 1) # thanks DEUCE1957, see issue https://github.com/Grid2op/chronix2grid/issues/83 stacked_solar_pattern = solar_pattern for i in range(N_repet - 1): stacked_solar_pattern = np.append(stacked_solar_pattern, solar_pattern) From 7fd473df51e433bd28befc3362b3ab7d6e12fca0 Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Fri, 17 Oct 2025 10:31:24 +0200 Subject: [PATCH 8/8] change the mode to be able to synch chronix2grid with recent python Signed-off-by: DONNOT Benjamin --- getting_started/experiments/quizz_wind_dist.png | Bin 0 -> 50587 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100755 getting_started/experiments/quizz_wind_dist.png diff --git a/getting_started/experiments/quizz_wind_dist.png b/getting_started/experiments/quizz_wind_dist.png new file mode 100755 index 0000000000000000000000000000000000000000..ad126abc2d5094e73431ccc86b804402932a4995 GIT binary patch literal 50587 zcmeFZWmMJO+BXVFH%RxQr9+ewmLMhFC7ptFH%KX+5=tYnXpnADT0k0<2I-RSGr9NP zuinpo-u>l_aXy@}KZ)r2ueIj9<`ut~ZTZ2jXAlBpaCIk)dCWq#{t9#?|K*zaxxNCc`ZhR5a$M66TwkAz-c4IDDC zJlqY$yPnBX{#v;>r!MCgg5+n^IizS(;Bv}sZYVWYq|ukmn0zpa{A$(*vuN!)>N8qe zgjHw!zp~L?&NJLZARb4#S*MvHj%`^0GUbd(QM z44aG?0sDgS7=O#LZ#2G(=fYT&fc+cd$3GsOQa|Z@>!hh%q~h?|?G1SeL`lo%&9)DT zOBI>oCa*V1X5pLGV#63MY@`8kvx%p-jA|oa`zFDk+b{#ADrzU+$T~bM`dMY z&6zo%|1*mxSNeYmhT`*I**(9-HKDz=@cS{4sMX?&!-%yNwn-Q!@p9nii3?iVrDE4c zno}oMz}3~2k+t>xrHgOpW$Rr>kmrwyzedj7@F@87>hSm|)oEgH>C3SrhA9xK;sjjX zLoh^Q)za6P3RyA#Siq~jcAqP2Ih-sJmAzn?%nbE?S?jehHj zpTq>&<=IXjgYn_VbDclS%^N)U5U7FHYPi|Xqd?s5LG1TLks5y1{Nome*~=J8|4rf^X;sEMa{6q{Z;AI?|3PcMgS&^))_8}6{UX@j{jj(*>~ zq+{n&F{Z5--lgbW2d5}s*uIf->nX`7QE_Jq#^bTpY}tKARW8AYWpCTs;8seeAkQvs zZMK37+TsjC)*e8VY6b?G4YJ|PIj4j2hxBRW`*PR6wZOv?Ibp54IZWxv!gg2+4SK<4 zds3!yN$5MN;r8ow(GASd;>a7c$&6zssVCyA5zBMS6t9i`$$qmymxerbqhDj zHdDTTIyCsuVt-FaGeUE2Vujfr!}(4ot4jG+GRLRw z(GLzM48L?W^Vj(bepu?pIHu?wE6LN7G@O_+P%3@qdwjWn0U1(%cDk^)XE#=rcO1`G zd|DTg!#!m9b>d!PiuXGDmdkJR0S?8moLK5;Wbihxn6Y_H|-aD6LnX^QP>N0SWR7 zi(^xIduM->EBAtNPnzwesY048iR&ZggdAQ=rARUb_;@s$+z^7_HsE>xl7)TbnlOK%roA#J8#&#vaE8xj0lpWfUtID_K~# zKC!mnYzTPx^VKHNZCT57VYZ5do-Vd}DAg>)-Iu2=XNB|h>QAa%Z%|#L@sXpFD-oZ} zjB5^!#9QctqfTl2!`GkDr7x5mA+5!plG0NV`jIF2XmO!fg4|(cjGZ*vS>wmH zNiOY@#`ZX2dAx-&`C5S66WT+g^!Ut?(qlK;ChemgvA!sWexpr7_hdSs$4(rmBgfy} z7JbuY%QfjfptqOLE%VHE2?zcXZaq8*4(j?fS6&+(_-*tZ>~wB!JU;ew&)Rt{+HL1Y zYqw*gkBH5ViVxzMH1;}K%)$%2;@jN6&n-{YnxjA3xgZwrV*l*F`&E76WO3PoaJN@& ze5b#wa0sgBPMaaG-bbjcs%0ij?_N<+(-pP4tZtTLg!<5KTBABoA^PN4xhk7A-J_+W zpimYMm*mYch7@g7noaUVmpe%a6a+nK2=teSPSIDxG1y92jeJp&#ZE?GzFJ$Oy`v}r zXM*g!D4LTWCAaRhdH>3*-P`v=H>$teT&+4d@k~P;|K0xcklRQV1Mfh;_9Qsd^|RZ; z%%?558`!n|wQ4kwfWRPV{Lf-+Zc4e9G2uf`!lizDEe!Fpd0(v(A!YhooD1toD_&_m zx9Ny>Np*x?;*|0W+$#}NN~#np6Os^VK7YigtQNO>0pfHEklca z_25~BjQY{$FxNnBl-mA)UIS^PhJu>Wo85I$PY!xDff}=S7b8K?v?q*4Ngu%9emrm z9sNugd#W=Pe4$D#XJ>;ri_Xj`Y4k)SEQlPz<;sP2_lMtWy&}fV)@QO>jN^id{X|`R zKE7k2P=hUT(^7^+mP5uX*q(K~LZ><|)874Ah&7$em!Q(y+Px=F(`o{_jVy`BidSF? zLT($3I;Rxj>Q26+njEB2Py4BIj4Z5+<}ez>iC>YPbS8VnSc86Kw&4=><_3OIB>Wef zkbq#&`6Ycj~s`Klw-*J>%JMAtqHHp&+dQZjQoNf+Y zGktSgq2R055Nlq540kJLC^zUUsxLk?b9ZZ=2EjkJB}Gs`x{RlRdQ?5c(uph0Lh!MJ zO+AeM4e_M@yp)Gszo6Q6Bj4z4nwo6iwRB=4;`!f(UBf{H%z)RCb0aQ|+`d1b6XHv& zOYHAI%aeV$#iJts*6Ky$hl?a*w6NElL}dDwklWN&%T2gz)?h*en;92^wNp0Cin*|D z8rQQGB*G6Qs%WNis7I6YLfjs^n@UFSC9x$1KF`p6oqOmTB}Uz3=b;}ysiZ@L@u+St5J9%&|sL0*cTNf z^>X~HnCHYEA~kGOMK}R16QMmXWLe%hkK(Voe9UkuN>q9JdN+m(CLK9wIh|Q4@;R$1G+5|Pspu2O?GyZ z1uH{6qfB?Q8LV&HTj(s*m9EdLGE#rC%r)qJa$bC!)x5pQG<36|Uou?z^dl}L#`gLA zo0gs70*xNZFQ&6aI*GFD;t#~#OoTho16oeB2`8M(pM>7^6oGZO!}suC$b3xw;kCoU zAijf|Ql^6P$xh>*s2irx)x@u>a`r!qx3MZH-dr~4*&Kx2Ol||pw}P+g-bKk}=`?u~ zn_gO-SIP0lxtQxbE-voK+!w>>WdGWU9)L|D6ir8di^+X5LVfPITw`JMDHy5(p-V){ zcO6g{iroKNT z0iI?FSb(!SO_w+K;|bN%kJ_3>G(P-(UnpIXG3@)*2+o&S6Gf(ZL(?Z9>Y|m%G2Z7K zF+%TC(HT@72lboY*7w&Gz{%!S_j-z|a5G40L1K))_;Vn}6)NLX^b+^E$`|F>P3A#R zG1PZ*5xUAJ5JyVv=v3`ZX=*fBKCe~1Cz;< z2&@28c`iicknfy@d1!WS1?5Dnk55Q3i)_>~*51uZ6>Yq)EJ_WMPSc%~tUwSLztlmH zF8qA%F@<^jp~uUhV}Zfit>~pCK6>wb?^>!F1Qe(Sdo7fsrE%RJivZz!zt|AV5Cwft z&~u;5vb-p;uI*AnS_LCFxd zf(_0xLyTPEeH5xcCm1e^y9q{8kw_?!C5XDxyzT1kQ&j^@@yvyL?1{wIwNycB6FnD< zzg8!whG~m&rE3b4%9+j}qwEUljWC#k^X)I&I<8N}jgT$;-j3leHH3`|=0w5veCyt& zg5&Qv>WfjuviHYv3%<7%#uW1L*l?Ux75W6n3JG}995R}NPe@6m_1KTix|YhpFX6*_*O%aY}s$Q_Fn5XT@>AeJD>*B^2w z_8{@0$~68r=W^Qd#>@z7{pxoyix2}tGo5U1NX*gG{%U%tm2#``oAT!@570HY2(s~| zU#IrTU&0J!JKh@jX1o>O$}g&}>=0UclKj|$NbYK<)hd5~w>Q~~>k%SLp^4(5o{;m?#2^>^O(dKE({GL`JX2-XywwX$i;K8J`p2xrog~rg zr~wz48NuO`y9F_tZas%OpCmQa*bNdAp+B!~t)|?5Ps#zPNF)gGRz5`9ds0MldJkP_ z(AWC2hF|-&w>mC!%(Wb%RTvh0pfk`h+rc;|Hr8vyZ|Y+#=Hy!}LHlGrx7U+M7w_A) z(*IsfsmGa2yEW%o+OAdY$DhZyMY%D4U*GwSKLRD@6msUXJ~8wSNwa}?RDM#Y)-Xut z>M;x?r>9zIiX$IovRa9!1)jTZ#zAFSlxXiG_7I^OC>v(WVcUZ`Q&}8qR_ki6S_&SbN}{-8@}tuaMl54o-jWemhGfoX@+DXV z-Wbj)VT_{~bzbn`u7TvgN++n&@KuT;18oWU~9T)Ow&Q&qbc5(@OOJ>?wv> z?fN5~U?=>BP1mXmpCwxLCye6`JcP1&EPZ4727 zs!52pw0GX09W=Ptz$f($fn$?mdUY$Bl$GnCw$w^!4NE8r>benYMTw?jdTDARM)R_2 zZJ5;b>q>FWq(b_`*oKL#Rb6ghX^U##2#z<>aZ_kj-6tW<&1o~UQbd&RkQ7|*z#pz+ z%Q((2ZdAg`U`Gr}B!zMz$kX!ljqlaVlrW>yMY@Y6-m6t>29lX|=#g9u1e06N zdkK?`u>}$xv(V}bd)+~GjDzR8J8%pJ>ycOipI*qi#>Bg}6nvaGLfB{$vb;)VjGsM4)MMG6Pw`XsdFv77ZGqvRb|dl8j4~peb#-QN5v7?U#-_K=^K%I z6IirAPd{k6XQjyzCATqJu&{A7Sbqm24v2ukYHrU*ow@3?NBg*J5#L>nU6-U;@X^>_ z0!aVTf~;9>D<0M~Lb1IZ!b|k<68qJL5ruP8RHx`gn`dnpJ+b zN|idTIqpb9&kKB4ZB5?DA=P5x5sZG>4Lg7>D2Sm1Goy@HvrMez2Qj*N?ytCLnpN z20jLWPQ_3rcUEY{WZa=Cdg1mDLT?lfJ-+BB6aXE=TF1R`o^W}O7{6$6kn%~FMv9Ol zf6s5hH(@JOFkyaCD8a5~^G#7c+Ocs7TYF9BHU4>VLd*kzyKH)fKw_?DVe*$_epJG~ zQFRX6YKb{=Ib-TNA!fD_QqDop%9|uMHQ)U+Bt&GyOw(_|BqK$H3Uzv&Z+gRHX|;{3 zk&`JZx6;(Kud5uh-*pGsJotkoGsH(AB1A|D;um@f&r2L5kzmk#g3Udyz|!`7a_?i_ zP_L&{gShK>yhPZ?&0>W_e977RV%F8;O?FdHdK~m;p>6RNoJ#c@Wmc(4L3&+O_zC0Z z4qTEr+m=U%*4xz{RDX7vEx}++w#yov%=IV5kah`u^_0W9x2L1EQhz+d&!+8I_q4)h zgCC%h)G9u-s_d)6UYGn34P~XR$r@J1G~K2^?IwcB4Z>c7)vHW@9VsAQZEWYUn(ET0SyeC zwZ8OkQ|7|a2w8}x4JB+ZIE1$3spjde~|MsuYM!=5lgft=i zZ(oM%zhCtq8};8^{YTXOEr8iKKS}&c3-E72`QKpv$NBs3x&3F>`G1+BB107Q&#amd z-v5Ew-sQjYR{;-YmbD+j?WH4B(1YAd_fVWB(*Kph-!iOL6UcJmeu${JvLI7F6yrr4 zLGAhZw|waT1*G4~Gbw+#|9?HkyB~Kw#LF=Eir1`MjG|GRZ?TtvixiYFve4I%7)j`( z{J-gJ*z^~4L!gZ(NDHGD%b}S%64n4)cRG-}f8P)_*8xg`@DkC~>GyO4R=I5n)EA#i zDZ{m!_?|d!tV#<7;%r`h&UwSh09DKJa>q~3j+b-yd)^7y6ilM|9vo`^IM|GQ+?T|R zyE=)v=lj~=PWDr~DVB^3ow z3pX#K5%;%~RVo`{3f@`dbl(wLkC`_P*H;?PTI;>pV<^L3_92 zkm?;xKgU(y7FaVnhvA$yMpuAF-1(LA{hFFbb@pEo?9XRWrTcUb&yo2t8vMo(_NPuzYe4PQzXyp5kNJN=UU! z50okUW{V{tj!~iDBAHu{McC&AVb4e84tL3wjYnu_4G{6Z3i6Q}TpEAx1{*qeG?UOz z8pxXCo!j#+zZ@n7B84|AGWQp&GBRNMv3+>}HFN+%32!y{1irqypW6jy{LVMoDVDp%0SmY|Qr+l&eL@+watkj|iKuRB^WE#kwSaD3k0xvI zaeD`@^l~@*T{18qzFM5i#QEf$h%k^@yZtvE*nnaiJ_rGq2|LkZv83XZ5m=_#YFD$R zQqCtBV=PVk^0}Snl7pTk&BrLYS+9tfC{COyffz$-)&{V2M z?xi~F$~D7>%9={_3%$8*>9`&lYJ^#zGXkbTeTmLgGeCdj+=$q0u4IX&ytydDxI@X_ zjK9apH2=6Ho$)9*^_|PcjLpLWWyWWy&PPXP(#J~YaY!^zhx#Hf$HqRLVLONmuQ%eP zNg{Ng*|R68B3W+LP^#I`BjV>Y2~;tYquJkyK|p{no?J84c^zTm5eIJF^4#$sooxSmI&G6y< z@*v#xRckax?8W?d)ezoxK#dl_WUPEr!sruOF*J9fNfXr$=UsZoUs@%Djjgyi^&9UE z5tk`fs6AyfXhi|klg2o!BZZKY1QL%z6gJCQEDJFjFMUuBW1$Ga0K4=LPloHES>y6KZroY;BV;a=^l0W3zlfoV-f;i(hhV}M7|K?EabH+Ok z5J_cUTr6|{f?61HrhNGgOe0SG_4}Pmg#s=HcZL|XN$D>UC6SSWim$8C88B!8`vcK; z`nH)?pjVGb{~}Htp%yUtYN+&t>W?NSU`xH-r}y5K&PHe;K2 zUqL#g`$x!sfW(ku;{Un0h7$9qD7H8`kR3!%SzA+1S`dOUcIBFy?VG z?6(KZs@E&kkeYo?^@DMO4}U7z@B$eXb!DDr*@o*hr9P)Y4s`d`1+Fy02a9x0_2V## zO514h7hVrb^*BBs5;wJmT;ovPeKZR2(O@&rrSmQY+)?QlCm+;W!}jo(n%s|(0kaL= zB#;qr`uXFdE?`p_2t#@+eD)knFL@~CU=bV!t@pHFY3EJ!suh=tPGAut1iPg#zPG;l zRRSZE$=mIPOvd-6hbpGiNmPDj)}$51=MNF3ldw=<&4%t4_bhyPj^u`a z9gf0a}YUK%$U&kz|#z zaF0a!4J~|t9O_FB0sx8+0m6+=1}%n|51gu#^!8!vZw3LY=tuu1&bJ40s#W3Kt^_C& z5Ow0K{fS5~RT(sGQpX;I0#VuiT(_n=0w5glYl6TsG!N?55SEZ%6C-T+azaKoFeVR{ zR>}D!INi1A7e{}X>3*w`;UoOTjmR?^%@zVR7qpLEe zSD8?l=ep1mo|r#hq15-?OErAV-lL*S%1xQ=W{0IL0I>;FSihup!7&5%GW+>7trb?gI+ZJ^Xio8C6rK=+%*`%!~a_=%)UXE>UQWV@?d@RW$|`_qC0@>cA`je4UQ2h#@m-hrg+w*%R7E1_WXsiJb zS$IGvX^H%m6d|RZrXu{0j?$dn@GT?;GuhUye8~64syTh6)I&X3_8BB>hH`sFqB`eY zaZ)@)k0pU6PD&mp9!|`4#OnK})KnO~YJ)<(ouG@0(H)`5HM#St@2#UU;`u(AIc!NQ1F_lQzA8o7qP+nnqyGE z0nJNDsLpgKlgPNmcPjwh)>M;y7Sd?7Bj79L8;Q- z_Jg3G{Ioi&0VWV`i5WsO(cg9mDZUo$YWNY3=vf#WTrw17a)t6WY-o zwrR(Mn50yoc}+4Oqu8u$!a~)a_Ik=zrmF&chkHUok5D80f^a;0>YD?AU<&S<%1qLS z7O1AkfuA}409=x*+iEcb2FP}Wwy%Je&K_v)s*xHOAMNKkH}X^~bigjgWw+d{!Uu0; zhUvG+64_G83tQyu!lSmp!113R3#OAQDiPPY2t?3eQ-Bjl2LCm)-UxNYp+o0bP$COX zs|s@rQ!?gnS;qx5@p&R*+w-TiVD4F#{%0s+iD$lNwOMd!M=#B}2FOqQ+27(akFyQ8 z=cxP5lPi(G-c0HrZzdp4!%IT%8v7e&Yb(CPjgI##4$bf8s%1m&|9hx|M;3=hr#+^- zn=4Sl|54`uxBU&UBSkpZa|lmJv#I$2rMQ%?z8gmf?_MVL*Eqs!`G0$v&yPy{0BTW^ z9pYMB*$}9t^}#dMzlM@$cSFg;z@LCX^v=*-TokdMV#0hDLI!VcANu!ZuchP;y0c}WcZT)u_xV4e)C=SIcE3`$H?uj9On zma|KgOvCh7v767HL56Q~Y^7r#;yK>NJW!d8*zRq_{vej;To^-VebMu%dg_oAgo5XP z*c~xm1PO+GVqzkuoWY~= znfu8V!(8*J1Va8%wWn!+g2t@heg<*r z4{RV(gN%!)6so~2seX8^1SFL;jq_Qlb6^YU>I$EF1~*F$W&?Sl`Sn$|`5IWJG0IJ)lJ&5;n*--VobNx$+wIO<<2zlz$xWI#TP%@JosW z;kw!i1%0JA1fjlaar3JezvC}_Abs}_QaGqM&&cUXd6FX{tW#;8R4MbJr?+_iUd6Zv zh@6jLZSGD6$XSLlQxboa@+zXcHXBG%VJwnjWCCCii-=zm+tI!A ziybJ|p9ez=RT(8(K4_2at*k|Cs1!cW1VdXPj6@)p$6~RMa($%X5yyO8xED@F)=APl z=rG249Vvu;sT2+h5_LzS~Xc-p*9`d1lWubkqqC z&KkghE+dz4NB&PuFD^o_b6M@)T#ol!1gqj)E2%sl8PrccduZxUOe}~AMi3E^J0PV` z$R;`R#WIctAH?$NA|DIQ4rPe=bUmou0!&zx@3z@f^As|d%NrAy>pnB|qWN>6%Qvy$ z5QoJU8-Bcq?PlObU;Bzjix^fw?UZEuv>!R~2;+&}ycflEM>2Me&7K3$v|>;52>5yI zM8jA%2Cwtq-0#*<=l&=ge^rzSb=<2gIJ%F+$(JVQ`)(PGn9|8KB49CKIKe9xg=9+y zciGh<9zJ{KklS zEh2T*CM?-BcD!K+TRo6ovLJM!;*g-S%aRDW@~Gna4!|l$TL-UghPMo!EtJec5>) zG7PxlnaLiLpc<|BJ^yR^78VrCWV+QU5|USCxu}Vl)Zp_SNCZ0PK&A387sOf}iG?ed>S^1yf` zOF1!QlQ`ZH1FG!>9T>Y{o{|%*XfxBS>|w-xg+viGH5JN!Q3cc^j1FQ<)r`ZuWtvlU zhi}jmF>BwI%@;s5N6`HaT_4`~=CGm2y$qL>ahzuVe%eihnk7eWL$ydbEZj@pPe)Y-}vKeoc z21yYi(ihcl^$C3vI@ksVO%+DP>e9q{s4Nawgxgv&7u1C}k*HE{_R7uG#_2l>r1ckuqR^!~32 zCIOe33)d=?Q^2(3y@^F>rgBN8`4iWWG{3E#KrA7g9G4tW60lS8@Ki&C@EX4VN zDeaw|e|CK&40rFCBS<38)T<-fP^erP=?M@P+{v5P+2-sGvzBu+Q>&z6QnjkMBHzaq z>$?ftk6d>qe2&S4BkowcCPI*$G=9zqGm2%dquOULll>HA+FnjMIMhORl6OcpHqn1bd|+2R&!#+hfj2VOq5MB=75I=7o!x`) zV?i4e8AzCjTW}FyvZ?NtS6Ff|J}vwRG>TtBcnkvnXb$5*1Ze+SKb$iGy2BXB9nc98 zQrW$>%cPDHdOUKr`HG4hd_H6Gl(MJ4k$A&ign1rBvV#@5-1Wx%o=DRg2fFbPCe=2z zzgkhW|IqBTn4u`*6(+@fryITLT!;^|M7^qfYocaxVrFjmcMk%SImYOv{{= z;D1rGcZ3XX-dzWHYj{NQiwTy%8Znb78?-*)a*A?hHYK;MBLLA6oVHm|Q!CSA8{0x` zJ!ayV^Qjo6+P=_@0hyGCryY;6?631G#s3b@f~4%YCwK z0dHh-x0~>3zbky3NlafzI7n16`RrPu5=Cx{J$a(B{8{n8&Q8Y|oGJN{Y98HU>9^W` z9nSV;5cd8>5!l7i#&>!N+zlz7xQ8JB&u}A##jf!@Cr-revB}07NiaesV1PdQBZp+R zJ;pEpgxe!Z0|zjrd1;~Tq{tH7p35|xrOy?CJ`|1id$*mTiu;GBp}?S_{D_AUvB7mq zByZ3l^!&bG)Z2= z>hPUydsgJ=F=ROEnyCWFkmM-8Js~jMN;griw6t}8rGFnKI-`3^RMYSj_LgPG1eKt0 zja^bPDhijmZLgI>=R5LsN+Oh7@KqE4@9>91GDL$6kH0_VcD@ob?ca_73_!X243t99 z0e5XLq{N7t9<_S^$~2?;XPD zW%sy^yu)J2d71XBV7Lwk7}ER!ifIhk2@-x+5e3`#fprechP|nJ?%RvOPoL--cND8n z0gw+K9-djH4JqfMIxsR4Yo*aK53|tWXpF(V**h>!!qMc}ylb7$t<0O}nE;nI-rCLa zPb7SSvHIEH_N$S_>1zLIl7&F_4N4iJW63y{k04)&TC1@4(pzM{=Q*lHGpPD#M;ALn zm}ElKK!hb>2&I{MSQ@5*#Tx24*>inNE`Y+R-zrhMO}e;ziO0uU4-3obpWHKENnmh4 z+-*kB2)X$|2FL(0AeKLv_>SVC$FuCGdJ*cO4+LXb`qteo_osj}$2<}XHc{`fh9@bf z9ids)JbvuAJB7OXd?k8>zweksuQ6xW<9nsUx0?+^2S`jp<#i0-lV2K^)U=2?RjY8I zi6e*o*j_EFEtts`Re1IO;~a_O!`&mH3pcwc+H+q+RDzvo(=kpGAbi_%Vy%;}c{)=G z-=jol`J7ZdvmhEx_m5dwsf`344DDV9_ZahvPf4r&oAOINZ3ZcuBgYrtidU1AMOIn~?LSI}`%x!iqX)Y#)->Jk1_7NxG`NT#5 z@Z=O-dj32{SAT!n`L?90!s<}u+PNI9>Ne!G@swn|H4wehs+KfJkClm5qVgGps>DEQ zAjD~NWqr1=0^Katl`1Os`g;-3HFefnX(0hiUJxRG6|tQPeE+U^eueFKBOeAO_ccgEk<;yQ)=`sx#pTy7h%LX~e1#uii ztjYHytrE>pO^?=IY*;hGl}Gl&;i=!bJgC_PEWVQY%B0pPq&9GOM;GRgUS& zBGmpoW`*IdcXlhnbY$`t0Gk*2T;A_XaPuw3Kj36u&K-;BYot3vVt@C;buLUGlUHs} zrT8^NsXGzWeV7zW)eX#)h1N%87fb-!d*k^{m>~i(oWvJ!LM0N8p*j1l^regX!Km99 zg{iYbFlZtTFFIVbo}9W}*Wr-X;>Hcikf!cQ-Ael*6S*3oWcUx&r8axyey@bz9~(m1 zyAov$dv$_#Ouv_L4@C-*#z}&;?jfq@p~Q*%8axFVp+8$<$-a1Vru8-@iNCvyf82H8 z&+YU}R$Rp)nTJZ>|p+mhOTsYczbE=gzd{!GK_l%W9+!&`P!X=SB2$xuD359 zzM!iV5zhEs?fntBA-F>#7eo`0uWpU4-MG*K1drbb-P_75}n*NjD(mPen= zfV(<4mFxtlg#kxN9{;nn3mAccqgVK3|5-S4;_ycoN#e=Cx?NxD zo%iFzp-}t|3xZM`SjeyL5aUq*uWH%8Z_qjMsXC|VVD za$f+OE#9=0Unik&aUeax+tt2jKAl%%zoLat;hW|$1LzmxpJpV6^<{X3ciOC4Uu5~dr#^(mqU z+3sZZxbh2@;!RJkoc{q_6TCqS_#9<)wPhvOq(m|q@fEkgskL{i)fTq-W#UP9`Q-xX z0Ldt;S{W_6GOT(%8kpq^%oH$j&h^ToupBmk4SD?=LbV@7M2$SG&??aYXz>7L@i6Ye zOBF)iS4#%c@hmAKh%6sCg?wMB$m_XYbR4_yEr%oX^{tE32#8uUS^!NnL&vwb|KvKX zw5lRJoMu@P%5};rUvRKoK<%CNvkI%0G)*rac*$xTT0Hynh1J}7>7#Aei zRN=SKi}Qi)tMlH75)PkCl$$Z^B%XFLqN}anMCx1in+4 z=7wzeeJA$a8*MLTN57cFOLBcWZt?it_k}MEtaVWb;f^t;h>sPse$RZF&xl)eUg$!W;{%&aXF&sX>Zz;9}j3(l7{DDJp=QlR^xGIs<|1z*#klGOke>uyvc3Z zyrJih-^}Z6brT2Ul$dbR$n5~BN)cZGaEnGhYQ9l2N(wUD@g*(mcwCqs7hGpnyU*MB zb8#xAU#reoNMnEuHt`E>Pk2`C;Vg#6s}U5ZoAr{sKqO@3611^;b;@eoYBkd8JzlCt z4Y{_O=5FBnrzMl}pd=+vNy#l52rGRIVGU?eLy+n7SnYZuM?Os%5CF|stO5Y3=dtcJB&=)730ppPWYoo zP`Y6mfQm=de8lOT26P4r+MK=5cZM0@#}`D27LLxhq?6n5t{Wb87Jng!d~Z^uo>o=g zK|Dh#hyfhF#sCG`k_|n)d_6wWSjVvAx+{v;_0KfN-S5?;&j|(zS3&e^SdNh>xBw z=nm@wS{g_8JmjfvoebKG51ikpcJz29ninPI^}@uvISZVX<}iVOQSgWgoXyn{wEzY5%DsoIh^NOLGug2EHNN`l^+VS1kfYu9G#FU z2tK~p6{jg&n%=wADtzoVK=c5l^rk^w$mrgjGIl_E8oJI?QI7z)YSyG6sM{l*A{SR< zI>d&XWvnv@ZsBCU)Q6%2y6SWsz^gvt5Kj4BUfUuN!FX#bDM$2eyk54}aHGMG%j?yi z()YIa`2KPw;BpBw5~1LUzS8A? zsK44Y@6fQsH%v-2I(ImRsVVBwj#wtdtuMPToJ~4a2=}cX!s3=0QtjJE(yEMbq}CyN z)!Al7mfH7w@bd!!t#$@!|5J=;KzQ~J48-g=20I4(_bb)P1DZ{8Q}<2LlOyf`8u51A z_3t6j^V?qJq?#MdfLDAXJ&O-PiwRJg2==n)zwws+j61%v=H?eE^*pI4S#|z2VEEI+ zEmD*3o^!7Q_YDL$Kk%jZ{fa-n&R;^v3)kp2AdOi#{)L79vzMw+;ZMNC>U84r$(>6Z z_*dw_{W`G?w$GSEPxi%X15_(sbr}Ms0=3O#buI?@YWy>p;B(!5M@Hj6b%Wsm86$9t z_VNbYeOY&iv0k_KJ!=Eca9kv7r3liS) z9Pa!>ihFmQAvJ)E?wL>eE*J6-XW$rUTmuFVF{LVGDwV!B5-Zq)(R=Auh zxAmcy20I!A$yuIkI*lBhCFApb|4lP9a^JM~* zufNlQyN9$Y&57@B*08SK&3ScjwLYE&&p4nHd@VJ5r{<$Ydcg%V5c*JWf4@0&D5xY7 zYg3amaxRT5^gB*r$HrqvZ2czRA_OYx`Dnm2@~JTlk6LgJXJ0Y{1PbrOVSE4Qj#tEo zhG)lLX8O?iSP7*jskdNWuK&Oly3#u>_V?d1k&Q7&DwJ4Btsg7`4Vujyj^7Yz2ynFO zTzN570+y1N)!bg>E^a}IVIN^QfdrbU|}C@fAgA>1J?|_24KZeVLI@EZ=UY zq$wADrQzda(Pgx_Q!Uvh*%@6Fq-#|eae~W2BGjuw+W`|vA?!+}`m?Nk=A6#Z%q$q( zuCf3mUP?J)+;?|-6viwU=@;6seV&-+cRg{e zqjLP}oAWapu&DQamKE}&zk%Sk z43#WAnnr%vk4z&TnMv8D5M3GUvjnwwD=&9@F?vSyw}?i+o15?M`&r-UXUUB<>As0S zb6WHyMCd&#l-p3r>P&i|t+H6PzL%yYY?Fx|K)-J}S_iW_J+$M!yI@93SJwofILVI} zHIGlhZ8!^nFZpr8*0%Ze(qVrJX&{5-NC?fyThu#9ys2(;Z84sh4xnDJO-6TA0%g&@ zvAvJ_{S3OuE|EmDl7$=kA=U_3Hzx3o=0gPtYyx-pQjz9((P~D?ZLSr9XdcZUVnJom z&JRVeb)IAc7e7o4Qb&Q>6SM`H)T*(4Tv{*Mrlo)uew2>_%2j;utNkaz(D0_R!uLf# zW8e#IC?579^39*!V0*e#bq%bTDT(@;3m=|7Z4b2`9(f?dNRQR*ul0=6(_r1hx@%N_4GIT=ufswIX>uiAP?{Ppzl{eFn6vWia>JyfT5smwlGjojUd0T8$mCmloU8a= zrml7mqW_r0^8$gUnDLxFp!0EQPK({(q1TpS;<_4V?O|-1}xNTn(-~J&E6kI;m0sQYfx%g!l$DZ_;H?R=>c2fN1VRH`IqT^jLB<_DE{fAk5ae`P12T*~z{ z{|6|9qry#>_*4bWd7*&(C);#{#3(Xti};H>$t^p&UYSmzJ!?Q;60^VOPKe?8Ko$`3 z;X7WlH%1%VmuYZe2GC3QWw5gTKZLyrJe2+SH=Z%d2r+iq8I*nB87W)V?AdpbE&Fbg zE&Hyrj1m$d$`Z0u#Hd6fYf)q`BK$v>7Vh8k{GZ=zyuP>ldyl!!b=LQJpYyq9RJ=S2 z7_GFk=wYTnDYZS5f)@&UQjvnn+f(mJc)gf!_@P_91ubp_ z!~0P1wT;crU$sVW$4S-9B>SH2WkNIUdd(We+rq?{wv!%ONyG~9Yu}CIrk1FZU9=r$ zuU7EE8=iI2OUclK72pZ-nsK}BQoPE?gIiJmNi}hlV0y25+;=+X`H-~y13iul;0Y1n zQ5ODKjH?4Mo;O~-s!Ros~0obg{;9c9j`>6=1Wf0${ z5Dar)4Kki_WmplOWM7jtjx0(v5ki^tNMu)yXOVI5w+}fjrIVzkX4Yhr&A<78HBCrk znSbv7s0K4kz!R+P0v>pP8wPTV*t^?DoPVdd9^Woqh~$3Eic7E;_xF5H)Y-iaV9pPc z7wbnRKj#PmBA_HABXn`$@mVe&urylO0MCfq@%t_c&Z$1|Rs68K$pHqd2|t2%w!Vmz ziY~-uRENE((5Cz*u*PVa%$nTQnh>jw{U*OT_~j>_#zjM-1jcG%__naJ?k&=y5AS*0 zldxy>GwEU9SGQ>%PSnq+Z-c!H(dN1LM=CGNHUU8T-aZ#`egw8jyi`}M_9rg)b&2c6 zi5R&IYDdyQj{}pruJ*kjo5z>X995sgt=G2p>E%raAD;6;S@xwLKAqmI6EUW=2cVx^ zYdD@&rz1SYgo2y12;fK}dS-x$CgzVGej&RDjofYyo))RAEBjt;)hYf&S*!mloscs+ zh(nz$^86Wb8OF%r!1W(;qZ{-PqEm5s(r0V0XHZq`7@0>wEB{eSiGt{=0Eh(YM+AhqgQzXHK}7wm_Z1 z)u)Lba?E@#2S8BS@dyya&{9=`L9sGx~VBDtv_M%>H^^-1Z9Y6 zhnKW5j-m&wCAfPZJ@^u%3(%y{%icIWG+2*QooOdCB%}=1b&Nv$;o8>)9XJIdwpT@z zd>@?M@<L?hSdkGR4>sfWOTZYxnvKB-bnQ?u}?C!vo~-d=OhzK7Ut z8TC;5HI<1mryiy@<&>|8=+Lxe(g_QA7#5 zKshN0Uo&D4@w!ja$ZyO#L+DPqYMfwQ+sitzR7SqMDE5;2M=bsn==TBxnOfJ{HLb%s z5nzqb%JqbYwG%&a`Ks-W!!GInJTiri!81X0KLJp7 zXI*%M!!9HYm2h<4-Ds>lQh~zliNH~?KF>Mx?q*1JpC^I{-wv&lIdBN? zK0_7|6k;;aO?(b}S!*Hrvpw5}g~j(f(P(DN%U|-B6!BClw{vA46}l&$fA?cb{=0V4 zWu>;-97a$AO<6$mI`vERw4@Osdb?B+lR4#SYodcnzJW&NJyP0%^FuA>qoaBe1$xO8 zpKhFRM?vi<0IV^1PxOHjiLN~S;2^KPiPc#=EV?a2wtP{~jL?tdZlzQ=DIh-)-bjk3 z>eqgvIC>2)+hSG_+|#npQ3J(DLH_U9%cI8vo`c4cFwj1H4x5F?JUG0-jWp?QO}D%M z_Mw4)Bcd4-D!DYKx+}Z!NsTu!h=eQmbnE+fs zZ4gG%?iS`6v{|x9%K4_%7y+LzN~f!q$v$In5ICP5j%EU5u+2&7@-~?DBIb&FL+(z*ldLC08mfU2Ro^xb;1-dGSg-Wr@%r`LLM3b3*hXq0!pO2L_?S7!} zudjF2Ck+f$zhMcbPFU9ZkDN3>&d)Zxs#CSyR_%`>_28ILr`^+4F-1@%V@0E9LH&{zst0V3`0t+wQ zt|#ZfwZsW!J4{e5OHc5uet0p(frkD2i8>dK?9mbKToCn!rK6D|PU=Y25o--6XGXi3 z#8Y=pm8CGCbjZx2OK?k}U*s_AwSMkkw%#E~@aXYrs5JC-)t7SjO_;cQ3UAAQLg^HbGp*cGMsRrEoc2Jn z;Qjd#Bm@E-PvNvEBqzhdQ?t*9WT1zKQ1*TfD3J)APQ5!VxAMIR7l@FnMjwiC%GL?x z7anBadU7hiwo!X=;I3x;>ojp+6i(@fkW-3#bp3*|9!%4qttyS)AhwWVsnsLh6ua?L zl|({i&+|gWMeJy4ta|=wb(UBPWOUjwMDp*Sa10U)7shq0J|;;#hDz>8|MC!Cn$%JN z7tU43mbBS5=H!O=G_?vU*ee0~nI@0O&w_Ua&W!;?niCS23I6tJVX0LATNJJ_CxY$P zzKc8-#%c%?Gt*`%sNcw!sF>W#4~&zEe`4-WwJUsY)$JNfu&;NCXLQwTbrv9%WiTn% zU(Em&;m+Wxv*b?O2_j3N&+dqBe~Z#V2Qs)RMdCk70?+MyRX0RYA^ii{Z9IOj>-$x@ zJ@ZE;mFW-lIefzQv=s8c6WQ%^DeVvt%zL4T+>3}62fkU%uf-#aY8cAaNNs&m(EG+( z;Z2QnPcP?#6Ez*n@<)_%3^5NWW64+;iFkWH9H~0Ju87cAK%i2q(ABwN5w>|G zoMlXX>R9^Q(t0qlNt#;)YTM zdQ1D#+_}AC#=cNr>e^djEp8NM*BDnX`M3Xo2MWpHk?&kz7^5g@P-hvw@Ke>mJaW2{ zk@7>%pqM^FML_Jvw{aH$UxhUAtnUP{b{z&2-;To-ihCPYpkwzu_a8^qfgJU=q{Y)# za>!AI)mh$TqW5m;`p;cG*09SQ!<+ec7vS@i6{ISeo5A?JU&hdhm`aJO}Gy)7*PoO{OC)9IE=G|gF~ zRbKIk$Ta#FdZ~h3E(&ck-EGr17;iwXgO$%f?}ywCzeP1Eqx5Fu*+%>S>f~Y>@XMC| z%8hGQrUJPL1`v3TO8t{jmtze*-700>hOM4O7(7dlLlU}S&N9$ubIZHjD<7ZYuI1kO zte=TK*-tQOy+C?b@}ZmB`zZp4$}9JZ;Hi@P9+H2f>KKCdB=BF`7aYq<69{1RJ{>(Z zmV%pe4Gr+z5=`s_y^UYx>ri)hQziU!?*? z$5Na+;+WwnzyQLJZ>;%l97M%i9hg`!2ry&}lADeMP%GvfI#;~DxQU2etUqf6(}k(! z=u74tx|K#UaBBo-ZtWNc2we`FV)$LSa7GvIeg+!Yii_kwwyKLeh9YAoLOL-3OChE7 zU{|;>&5s#x3cdT59Jv6+mCqGo*S-w2IW>)4x9bK#}h3i+?cTjNb&|IHIv_}Mcp zp!gW4N;2g72DLceMv~?^$4*@OkHk1X@JZWwtQPrI2N_ija*s2;t);2xW$}qq(N|By z%YQhc$UwZ)-E2RMw-!M{9T=T*o_|J^HW12&^bBT04@Q;PLoH%jQ-&(b316xf(~SR5 zJy7Jqya_k$FhXIJKb&u8AExtSQ@(Vd@q(ioCyoFbx1NnPy0Qgn{Btr=RRM)F+C`G- zCf^vTBvZ$(g+#u)(b7QWaqSYDMZXq$_M=sH{$o796KRsWZfz!buZI z>7sAItco~~C^QfoWD_wKH=5a)J>Qs}8y*wV>g4 zhmZb4Q7X{g_H$#1m0;$#Z=ZD-b$%l!C}Z9G$*K&Vy7g-AKiS?KREpOq%`P^~SIDJ< zr4RlNrLESMQ?^T(#C$T2!Z#iAqvt;vv<`6r2lgV6yP@!fF$Ws))@lGa6^bk?;!Lua zS&@My>n2Oambq+W3$eKC+m+nbPgW%2QbK}90*MR+623aMzz`OZZH7v|^q`}7!nN&d ziDgyTx}(RnN}h|8e&h<)<{Q3a{JJ)Mm7)GKk+C{GpX!o53D&bhqxD`~@o*i7#eLSg zpv|o9$pG>?3YXW)b^fV*{?E48q>`tv9j~sQCVv{n%bWC}VD~y~qaGP4p%;an^A+$8 zP>zjWw8gPaka`U#MD1x=CeW}>mtdJxzp3m*A{l^Dh?ZfZk@33v9SZ9kmo;m(;Ar3bw zG!itgz`M%}dtPoyDk4IF(rH*i={#tOjUz08pM73E+(2e55pDSDKBLS7-JCP+Vfox= z-9)qzh3;O0714*mu6hiLUwkNjIk_-eT=^=F+J4@raxN>v`1b4=M{ArDajUlH*Md&% zMT?yd>+C+6r-a)I*xD@?DA3{VcpkL&o*-L-hjvd6R49W~^!hc@n^atbQ+gaH8Pu^U zR0O0boteIyvL_FXm5jbefrUGJ&&n3WA_b)o`nTs^BDps*KGAM^utc3&;$un;WTnuv z+&;$qT;^n)yhU6cb7s06n-@(;LUXA3Y6rd&H>NZhO-`VOi)WeX8&<#er1yl$v$&KJD*NuQcyR0!&vHLkpFgrk?kWC zJ?AHw8(j|}7mf}>6$N%I$dZ&}C|qtHn!-V#btCqawa4=%w#xWkDSz*;rdxJXSIdl_ zAao~qx$%E~KJgIuz}8fD*N*I_d%FIy#-EqZa+;0Ax|`z*who?@dU@JIfiX%g@#9vP zPmF(A*cvD0!;2nmE;(CsSptk5FUKe(KQR*fd6~ZD>N=T)tUS9(ncCRhJ~Vb*$VhKx zqtkrrqbqlpl>fB)3!fTo;%*VG+?D&iR;!JnF9z$%4X3h)jaLW0Cs-ReVa&fVBW?Tm zMw6J=v~0cX2+O|NtLA%B>+vtOZb!^m`F2ptWtvBC))Vqd(j1#`&6^i7k0Q@J-5Dxh zQN|eg14Fj+bDme(pYUpqz2yF)$qloy2oglA4nwtKt$pWwFhBC-0Lg`6SH>v?y{Hxk zU5^)DEhiVQ%0j$}H9O7|zM_`zf8!GutD8CDRnwO>HwPc{=y)JT z_VXHS;-3^$x{z>a%SF(A*{Ku?m9^t48?+B#4b~H%Y^rXb+%*jvlQE6jnO&I57)FY; zMYW0$ZsFXS+13VFmbBvc1YDNBPuEBEKlstl zvqZ!AcMHcFqTC$st|+b%dXEi~hRpTg8Yvt-5B-A$&z8%|^@_t{ZO;xnUK0x%;jIGG z4Dyj;Y6=rvIjzx#M<^VsCWJ%sU{u*+re$G*6ptSfXuWS&+Gw;7*Y7`{I| z<;4nf&p(OY1l}P2!|1X38Xzoh~KM*YtV4M7RYr%FGAMj$`BGP(uDzZ@CmZiqGqFv4;_ac# zksW%g1)`0B=4n}ORqV$zr>jH*;EldD5fi*E?@a|ht6RX2CNzAcM~#qs%RI3rj6+8~ z966r|4|KmE$-sD3XyKO1dGY6q^6}aEKc?CVM1+DDc_b@8oix5IO0mC$QLC^WY<#e? z@cGpzt1q4&bNPh%O>>68o9KV#APPcZglI5Y3|L3PIWieA` zq+Nc;kIE2%oV-?fyt8+bLr3{Qg6%nDXi$thakn>QiY$QncdN0M7pPYdZOQAM^OXm?B<}iEa!9}Fu z$S3;a4sRP7uDa-!`Y-D67>&FT3u)M5J<^5>UIo}CDp$dxY3cw`^8A+R3vzzrGQlz}UK@dt=~o#rTFZw7z&I!j|kIzAWbV~OSJtsBgGB|ZsZd<SZc>g+J)fR3DlwY6^pg4{W2a~ik$aeZ^@XJHX$y+z%RyvOv zOPu&aPsNPT#5-`g;(VtJHhSZ`p1r`VydT6fFFvd5N;S|x_C<(3K#CS<7!h?;!{@;T zZmHf4xTC2BX_yGouyZOpe{($PJAotF!;_eMakR=Dn8Z=d@bWS`-kc$W(74Dvv>9lS zum9&BA!MEazg;X80DeI-gF6K;Y!QGh&!pkvuE!?m)FH>9J>nQ>JqRg%C0_Ye1c)KXwM2ftTGPu`JBqX&+1UgR#^N&F4vU^r|d)$I{Wk z#Hv_NpkcLH4PdzQh_{l+`5q#Iv8V_1=CIry6VizF4ss5hK$BAr6NTFKsByWw+MjL9n|vfaR? z@u@Gd+3^P}L*elAlpL)5Z17CqmiNgb*|}6nmAkgq0b``!XehA*8QhJxgyC|vK?S`c z>aups5!63@v!A(%njL?*VxTbI8AEaslR=Rzi$Zuk{=qIA%FC-Y(rRms_eQj7o7KacJOWvw zHQ~>H;voe;lg43k|IUd<=(~F}sKR`JarL>({FpTm9iOZ&^PPDEwFJcKSm}_N)0bhd zdZoB62F}d|S))WAzIn20bqM1GIf2)ArE7kSlLtYTn$39N^(?zrh|68Iu_CnU#c-LBJ5U>oar3_ywiMqW<^S<=ea`Y(dUb%Q(z=8mW+ZigG2! z;gEUGLHrM3fEW;>a>8V8?fS+aUGljxHlc)(_a{vZLGP8@#&grYeUC z4eKTL{q}ZxcTMTgE1??PgA?vr5YvZlT14eQ)WR^Nla(@riz>3grtasIDo9m`)L=ifRLtfkq_?>84cCqFWAirXIO<_^F; z(`?0UZOTzKI6HOhb&>x^7Da>meTzeKhvYL{dEPT6qlDP`0jyQ*dg#A6CCH7kDJ@cU z0iBhaZK%!yx#krCjQkQ~QXDEWd*?gtNDW@1rfh;y^(m{}rB%U;HBpsK&z{Nzbv_~` zwkl55lPi;M{1r^I*#6n}5v)Et^Swd~0#V_A*R#?xpus##8bu0X$_;=3h#5fUiAvkv1s9*S?>q6~@q(Ff@_@^1kx;Ve~$jb|~% zefc5p-~B!xlvK0c!;;zs2*sYd_1D%@OAf$H4_Q|3+;LR(;+-cgEF%4gU<=#Lg4bT zh~F7mSPM~D4klo!lCeiGdJGGnwnf7dp!Cn?_&58E(koR1f;Beth!OIDbc4m&5+z&* zM!5dR2oR1re}n_)1YC2^EHfU=34)JE#owI)luXrfR1ck`Wfsb2U<>A00+c4%oc>&z zPyB7DE}?Y#N^wN6;8> z8;sd>-c-%q(gd@H$hvr+JtSSTHZa zTemd+C*Ul3kQGq{nd|5sf+nHlgvo4KWe_VOaA(61=8qb|v zXJ=y2^dq@E4=;r~2_>GC;}$m%2c4wdUo*O*l6`_uH(%RJuGR;1g7L5K2?`F}uN);p z7qJ{5!&u*;mZL5Jm#JLqd?;pi5U0P(5SRg!!%Lr5Ky!^-*@xOChA@4BWAd;&l*u_9 zZvM>%`5ehQn(5v$%=c_|v63#4d&^m>IiLPhn}HAOw2Of6P(uR}jAK>mcW3~(B7xxG zJrg%7rT7@$kit&ksM`DaYv4e7LyYcnnitF8-aEV%+@x@}{^~iB9LTcp-zY&BM)wdd z{WAh)lx9Ewyk>9cd_xSY=k3pjA&CzzXLn3(8bwb4h2)jO1B3#hl+sYP`=AV^{%bT2 z?7*4o^wA;h01?3cYScCjjga7v6}iV-hMIkl35_!H4hPqmQ4mi<{V^)Z;rvjpDTIb7 z(97;B(#fawOk;)@t*#6wI>q?{xN*7GL+w8SmnqF9Ph!C{hfJuN@-MJUwqw(QH#Xb- z)1k%n}X8+vo7S!W@nJTMeClr zVJ>(T-z`om_GdUqeo+4fIRhwg%P1;Fn?7)*rDt{gVB)_a;kc>;W?Ha2#;D5q=HOSrYfUCz%we z3`r%!eF$P>;&O;io~CTwX*XK#U=$?1(N1m2&V_kN}JO2G!qC|dvOkw&- z$Yv&AV)TU%1J|L~BEdoPYAfc2YP7Igj^D)8mCt$Hgf{55~J~b0M3y zA+2;BQIRb#5vqBqH>6B<$e9$Vsks&>t16FyId(ca=nMvq@V5JP#_mkTJwFW^P{s!X z)%ce#CUw5?D6DPydP?rmKAFO+3nP`w==hPtNUcUaSdi zH6KIKB#nD6BLLlXoBkBXUk(Zfg@L(_9Wk;5L^F<2S!QOxD<$nkwy58eT zcl1gpeze^*xI^~+)pg<9?xB@LxzbNN-)J8%SPAS5vvO_Dys1XxXLJAJu`ZcFR+RvJ zrdaFj`T2wj^Qai1?Y!|DYzBzLpahEjU(TNu%lhOnX%?+cE^cw#9sXY~!3;m;wq2=M z^$@F7?0Z5f>*YyfSWnt%d_vpoOnqXa+En)92*R?Cn0nyD+*Jxi?vcN{0HTZQWkT-X zc2{VL$`-@W@l^pKgk@Ll>qJEi;>iO%c>9t@3c7bBCPU;I+>1#RH!T|{H_M1OO^xDL z4ff$>{W1x!Qdlb=XR(HCn`mfGMD*a#ScSf&9+$66?7u?pKWBEV!FBz{gg4!FUWvOG zni(xjWQdqBEk?WEH0=DcaB;w6=lfC@d9GB3$$d;$~iO;aU$j+;xISr;8o@2+&_bUkHA1qXH z*wae8x+lpGydblBiuY|jM~LStTsr=k{1=L(m0uhkU)Z`vbw-{X@7IuuxJ9=^My0^C zM)7XGv5{>2(!XiO)?&@O6z%Ac<_Afsg$3N;B*#k)k7ugWulx=PDHiH zbd~2ayhqEVH7Ectu<-oU?kxEV7+mam+0k;x-Eh}(xuHyapb-ShCFJ9ooHe4TxLYJP75A;wx7~sl zkUE^;oagz;D`3*n_zFx~mNnn&b7DpbAm*=wj3fCz3&z38<~##KoP)|D&h-n4cQaXH z?dvFPw}1ths`$4Q*$^wj%3mlXydr05j>VLt1w^muhGR+8V1-|kRgBKPA+o~3nB zvgT36I;O6S-QNNqycHZuh^sDMd13N}K@%Tp?4&0M!@Tc@jWRi*3 z;kP}A+We#%Ty+Q#-flF>H3TbRcy&fpm>{oFn`^UH{)u?hJ4@8hm;E|Ntdk&Mk((O^r!sW3yeC~@^N%b zhd%=>k{t)V&%D<5rPsn%F186lmnCoR_cHJ2wFU9CRXKDn(8~XyVi;K?;N`m=8@Aa} zI9%oQQGd&4{qr@!AMojGg8V0xzFeyekUCLiJ-l_?oK$?|D*mtPaY6fbGBvPyVs`0= zc}d;a;>yTO!c&H_$wp$Mz4X>L1XydN+c=p)qy3j-&(j)q&w2n3G=WBrHm{xbvpUqE94t?d!W*oQgLbk7JXxw2{Ja~q}SNB>tZDqfoC zWg+K(`PTcBf-RyBy6v#`dIk+gZn~I1OEB~RgzH3D+3O&oT{f^pdTs)xV-L7R(Iemw zr#^Y(KgiHBiFN_R1%srV6uE!{7mJOh;O*QR#Sf82g90M-%bX|L992$u=azOj-Es1t zmA2LDP^~*&QrV`!q#OAnX6Qov-*wQKEI{Me*D^tTxnXom5sAmu=_`%5Lyj%@B}Idg zpBeR>^>6y{&seoXFfL)`g=^%O6pY8|qP9QtmMZBO@^vL{uyteJVN6Bz%-8vs+!Dv< zIfirp5Tuxl$zl5sW?o7BHu`!{qYM@+vp?ZT+$F#B-w|`C*>#8I8Wft4X1<-iBAWKX z+U?%YXWJWa>5V8l2HBzMO+};pONiN~tuUHzytOz3|E_J^Uh|Cd~zIlo4 z-obFrrwlZPsavprJs}H57z9}`?O&$!suR-5A%)q(dy)0ycl;7=j}Z!|Kg+VvD*DQw z5qy;y>G4`L6+{6#S)0)FV9pZd95*9va0`r)06lyR=jTQ-cv^N*+W z7N5neM_0yc<)d`08_E`?BlqulOBpK}z!;k{*z+oasy`zq&i+d_MQR|~h_3fPa|?!& zz^&7_;iV5Ah=7R0wTWQn8>;bd+?Q1}99-#h#h8${pTQLIkXtcD=1IAdHH13dQ)el_ z?3^f8xjVzO8F-PU+!jcYnM-soy2*9WwFv3Uw9SV(yGiBjgn0Win4>}5>)h3a!T;Cw z^pVRt5EXvN z4qVQ%ku*zD(+I={ul=>6W@l*XggnriU4Qszq(jrjl$LEf11_^2{}&+f?n}w zs643RJ
xA3jpA|8jzYd$#`%UYKHmEQX_&Co3+aGah^fhW2`Mmk^@XGgzYCV0tCkc*FK{8r98*DP)1LAp=-O(qNzB z^vN&%9)_?lQSA4U56(V_@&mZ_qEcOSo$fR+(x{#g(K1mSG$PrD86NW`ig3IeETeGe_pX4XYnQ7UP!W<4y6 zW>n}!^$rDz^B|{6_2$>Exu}oqHF29qHOGk#>;43Dv(?VOtfeMk$4k^I*9$qJ>}~V< z5K+YI!5DA5A6Ib~YtR_f9}|s&BLkTk_my>L#MdFe82e)lY!|vH#^wggU4)={hxyJQ}>9! zm9>_uL*aQw5beW?+&1f4u1N=nxcs;`Um0ve`&r9)R6E`=7q^&kcGsVI0UdBzK;~T^ z5(@GTE+z{#IW4H*Bls)yyYk!%b=l84x|3xL9_V;-8YwDjapA(qcXC}u_y^0>bunbd zqj+GC-^$lew(Un!T$>sCEN7ph3zb@EV&%pBB2K*-BLAO2pHX z6?ZKiL_cbh-WW(E1xxR-gW`iLJw9DIF)7npDm=et$PKdNjX}kxfGrrdyikkvmy*Ft z30i!5)De z(o3Gv%fzVTcU;56o}MncN(Qg8@ z1e5|=D58j?)hIR$|6=+~8>v)-LPeP&j4}8Xv6lH; zoA0Q<^GIwG?;f znl^R~>6Gk@599Z1(Hq>I5>iI-S5e7BarNVG%7mC)dIlkFG>i?3`BSjJWoNeeTtmq; zFs}5*`jwbGv}qtSxue~4`=-gs2h>4!i2SsNly@0Cjgup%=u9hrb#VK4JZlB}^!JxD z;x;)pc98O){AeF)!3?M+#niALUwXthPn&}?UBQe<%g8@hN5fK~X`OPJK&A> zBmg(jC-SjFr`y5{pp!B+@vDI(84&dE6T1(8+~ADKHr$N!t0X&7^@xRhj4iZggz)>r zP(%#3>g!iOmaG9;14}LDu-1APdZRljvJ)ZW&&wMg=0FKu99-DI?W@^Y51X$zqml+% z%_fG$AJ~V}je%cWxg-N0vR;6V_taO11bMxBw^-E3pivHU`)n|)2r{w?;X$N?lS4w? zwOcdXxD5%lA^HR!==CkaNkk{o{R9&h{<5}x%qJavPGtNro7muwQthqAoYhb+#&9RonPzlSx~AXV?6?6>(bk2qa2uSxkk zXFfI9OW$j{C1y2qZBq!E$V5aOG_Nf^z!T{Lt_K7K9_-be{cGj}jVpqnET6vg=w-#g z`tui_o>>fX7CK-v=a=j6vCRKp3kz+osUE<$o;D8NU!KI2b9_^kEGOGVo?88k)m@44 zQ>ryl-%#R!Js!1f>r!-*xQzPIg5>e|C`iF$c!Sn`j_KXiWeZJW89y2tNkWng7FEewfGI}l85T{4HISg7M zoUc&o<-pw?Qy_10y%Ygjh}gyzI+q2LDH0P4SX0&}jg4ti5&W!*YBIXb&0}xbX_nzn z;`0F*P+!7A_BQ-0h`o~LS#C_nxH*&8-yz^#xO2=Wt8*wFJyopM9C*g0iPu$Uxtqum z(&yxW4LTMVoeDw_O$EGW%lO}H6*^<~`?bxExpChN%;>DYOK7_^F_c{sOOfY0;c^z5 zEuP$~QWzzDSXZbq^3mLwl=BHvHD>kw&^_bbH)e-swd1!D0J|U4`*)TW8+aBenDOw}y zh%XiVTC%#&cIAzLHsA&kLi|gaaqde`+aFR(pG$(iJ57aTLEIT|U$bWkMimNW5QT;? znoa%2#7l@9!@`yOfya66i2AP0bFLa4{~hyRph;-1*}B#p=*Zsry}wFGM;D#YOEnA_-hBi(!<)JWw(~`0@UyP`t`hT6{?EHVo8EhL(9rnHmb?bv@gX`v z{3>K8ApPU=WO0RSa>DEyF z^{RDBjDm~DU(`uwV~SW*rx9wEQpUyF8ywA0F?bt`*rU&0n(le@qGiSv1eq1s%ZpUssQX{;f)1 z8e0kkryy^P2io|M+iDV)nc{JAp{74A!p|?2M?*kvh_&(xJ_A??|FfHIo#{l9)fqsW zO&)=YZ6CZ+dm^~#j8LBqd+;?=>p(jeutZzgI4tiO^wpTLBjzV7Why};zfBoYq`8dP z9+Mw(-HVjfc@ECqf4QLZE(5%44IwYv8x?cp1Z(6zNp0)8NpLqtT zd!rdZRQvPH3EuV92xl820)#Jlb~9bC7C}M#$EWk+da-zGs>D54qFV(E^b)}-2eTU4tDSGj1KC(5hT@%zxMemRGo%2^@lkNMq z+NiYo*b>Zvkzu@IuRTEj3VZb?8I2Fhqcck-@kLk$)Qtyo(3`@(zkg@Kiz$Y`*5Uy5 z7pe9?QR0xek6TL~dQ=~lv9$Q)ioH`BSP6uwY=yB9PMdErSUy1OHV}uZxxi18&6#h!wO z7Lz&TnbW}tsK`Jpg}>C_n_TH&6P6o;^SwdT8@PDEO!uZ+slPGehVMNC?ho(5@nI%g zeL)n0Lw}VJHdlS`6!sP7olh#dyV-DZn~JW_;yFaf`swPC7Q>y@w`C;|s(=FK$Yo_u z4i^X4t0~ZL7y6!oZoBm8(+TGAAvxz@g?TbVOPR=aH|>^Jebfpx#Z&Kjf>y4?Df+*D zM^RsTXWlnd?gmuXOFyIv8v>TNqzGjX<_ZsP$AFFdn}yQ)tNIqChWYB)E3clmqsHly zRK6>Lvdx{qx-m^&>OajabMr!d7+c^64sL#vmapZ#{tCI9HH|H|u%ZV#0 ztBvdt933M;a*SR!(c@I@64fMQDi59G2X*8yydQfy9{&o%4d0~#tefsP(^~K7s#1C- zJy5k3ZQ%f-#k0gV5~P7VWB2;Y+$)I#{96f*(StM4>F7m7N8a!M*lsnKr>wvX@Y}Gtbv;%FsNuMDAf=KKCW5 zIJTHS_@E{ltSz2w)x-w!nXv3PUoPIhZIhkQccb+dMc&Ye8pH};@Q(%749h9!@LhE z9y7?-qjKdUh??EjuqxhV=8j*I^@a%VwG#x=vSn=-8N)giF1$BxTKhuuq0eYfE%m8p z;zrZ8a$njV(S5?qx&65O53)L7Z|BKn3~5p+t|4E2x&k{IeG)kz+h@6&Sgla>Vdt$p zHj2G}msdr(k)i4H*)tf96w}Xbihp+jX04DoZRlQa(O6`x9vnq@6P1p3IVBtWneK;2 z(2{hZ!-Z>7#&W*2Gti9#9Y`+HZ1+mYsZs*fpeT$@G~s z@V5hw??x$la42Y>slgyo9Dm434tuv{@czNB+RFtp`#5ptA=k;CVxo~cJcWck(?y>j z_Pwn0AIqdbfU)-2s)obiK(Juazm&0AE_|E+8{U@582*LY(|20{WQt8PB}NW#UVUN~ z8npQ}GXjC^*)G9v{CEa5LURyd^7{o)V^`NViULV@L^Od3ZC({YuXg-(mB!P~W7cIU zK`^>fI5P7Z;D_F5+pI$9*1_Uco-z4va}kD+P!_haNRHcz#prb`1%sIRC~vN9A0}l?;9QXEET^S2AS*dMDbz$;Z8YQ#Ey~I;~x;1v-mE=Mxy<6$lzFRU2Ud zD@#9oZ1M)@j&M$O#D-J9+e{dkE%1E;&%%Ow)*lz^>ttwUyGnG1s8<`HwPCL02Wr?M z^ZgqdsP^RjXGU`wXwb!L^X*;tx{2T5z8Gxe^@Uo_ITQV5XLkPkO=N!dKEYgd;yly% zgEFC-%JhXct@e1p6wktCoH%zGTqb`jp|6^2JCh;V6fggo5YNYh;1q zE|{QSLpFMIje5dnGF%b&MR29=7k;jJ?uFGA7*6AzPI`>84|-JG?x%gCMsn=ZC0nP0 z_fJ=@_*Bqmz8V6U?#mM*@8BF{$y3?&Wi@O8P=J`Ibr%;N3B+h<14hH|#LSY_7-Qc% zg9P>&ZEK5u#j7J?OgA}PYAjOYc3#_7iF$&Qm$ebg-Y8h#J)B4}$+_WEpP)Y>R-`js z9gH3HNbaD*MO9mp8g_{WB*!CAK;<-qq#a4}ZUqRpU&|w*gOm%urTSkxZ)u z9$avNFHSedE^0%}P;80RFb)3onipS}qi=UaRmfL9j5Jk6Yhhm=EyeJ&|F6C8@M>yn z`bG>=6e$Ww7Xj%Yz4s1M6A03a^s4kwA_$^LQ6V55q}M3DgY+Vu&>?h?Dm8$V_wcyh z=id8$|H1d&WUaH3taWnsnmsdn_RQ}$I~YXz`xL_PNn(zb^*4eEk9Gk#1}Wl}!gugl zz8ZGBQ^K7xIro{sl`pJ8jhAZ+1lY^rJ`X8TT> z_8R4zPIc{m1pe3k{2@#IXS>q!2yJpwBSOe4HenHA)J8oIw*{vm6#OO$yc?Ykt^lX{~cPdgmHJ|ZR;TefN8?{pnXo` z$U*f$ys0gg@T_Osf!`CpIG&tM?|lM_Cpb1_mNd^YnztWOp4DfVWIDZ*<^l%z^m^OU zos%j&*L58OVgHq#1~8c=;kGe%y&k;36uV!$z6Q|c{#@mlKM0Z9rA|tfYq~3$duWxz zRYD2$BpMw&5`7W=6$_1LAFC|u%z-49z(ky%?_HN9pd zz)f-*TYj;2ypA-F?b^kA6l&iBa+u!|-_irDnpt3DHB6|i#1J%fy(d=6i);G(*Z{0* zARE+5?IGG|^H-ZSy#%SdD}r?W3>@EHbsnkL{{aftfLPt3Rl6o^JY4PU7K@QtJ>YBi ztGh^|bzzwRoWVvd%~NXpLv)|heC3_?_eOd^rY~p39Dx%xcu&Rm-j}!cq|j3rMk_)G zEAg-GIf4jy>@24$QC*3{E9R1`i1`R-q1aq5H~!)LqHFI_1a#yng9*$p28O&r==vAn@;zuSvoWYSViYapO9h~qcPXa6%3f>KvCY=^`HX!&D&i`(t|65+U z=EikB0c4F*c*7o~hSMK^B^BZq*p&O7X>j7lwupXgm>ya)mNV#Xla)S0hfq9jFO@uW zhvOP9TpwFl58CH8cx}sb^!;t0WW^^LVbeyeg}|oM-y{N8nW0Shz=~StXtX;X>#dru z%zY6I5T_PH{M=FkvsY>28Gr^t11x&gxmp)qV!fh|Y!F@7r2hUca3zbb?91N#Bjw;up>%~e6jDu1JspEUI4_t53cs-C zb6U1x^?7_zvN~L$n%;N#qC#~R-E5Rb-U{2LqY*?Vl|JkT%#MC$dR_icP; zb{ZA`ESwxJ&{0XpbqUS;CFXy?XBth5<5)xKZ)+F#5NIh#ci>_*ylu9?S=CnH_a+<; zO#z8rM#`Q49)5s;r^Lm}@XxTDhE+;?15C5)T*+l)La22gZgcO9Y+*0ynJ&WCa5gLLfqz>rQz*yK#QY{5B)vAY8|w94&*J@n(cJhe+u`UzAn^H!YZ|%-)Q5&2 zbyDz^m|qU>y0MEAmJ}PA56$GS|6F2=`U~p(M@W1Se29KKxE;txkgi{q+-YTX^fo}h z2GoXU`uIK^S8{w2+6}Zm@=y=milaH+twMSW>a27O27}IkxXn^n8vHn@2fDXzh>ZYP zxhmjwK+8&zWB=~U&;B0FM1PQ!vwA2OTsT1I60o{f7SzVC|E)evoZ-`O}S#c@Scbu+_g5WA9{bh~n-`Ua+tU+iqXW z2!4!~1CV*1QuOaSoy2sp?e|=T{*(Wp1(cObL?f~6{>|MZMY~@Qwj>Yf;{wTp&|^Gd z2?@cZU4I}bfB2|?)qQ5x)&N!{he!E9Mfi3yxg_eu{t#`cm0o*D24Lh5SGsE z3N(z8NYV=cny2Lt_N02=9Wb5@1zNXZ-P@|~y<5!I0VM0^>l|5LG;i`H9+C%C(9;b+udJH<)jY=T@hMB1$t&No;01K%HOPVZuN39=rw#Ny z6*GW?&>jQ!zrBm}ukYbnTiDD#jBfi8jebnPyQ*zO3LQ}^K~(re!Jer{Q;ZvJcN zwzn-JC(p{;#<{}YjsqJ4c6uO<#vZvhA!;^*$#hKs4ph{5r?o^EmRUtWawCjG<)Q3+ zJy73#Pa15FL4pyacEgfJ8Jr)ME~>+NEzgaU-wcdPvxR5r;vcr|MrXE*vUq%caWj)Z zc&R(q<-xBqtpt@?PePGnM8X3`&UtSg>~`m50{756P@9}DLiqQA+WTpNk;vAuLMXib z`~XhSb)VR(KIS(rIZK4xfqT@tR2^uwlI-d`V_UNzKj-W@ak5zcS6COrxxW96p^w5UKY>CI5xI|3khikN?}U)Fs07;cKiqu?Ux6 z{x{5})uL7$=s(FK9`o#PvxMW zr}Ns`>15{>tR>k~;S~e3kMe6DAeaYB>FL70Qk&>qTMr2+>S0jrkES4veI}7)w~%QY zZ=Gwy>=2x9Kt(&ywhBwr@TBa*<3VSelKi}Y$*oV95rwN?6xE_bTwOq>n$)k%6>SB| zz*cg-TwOYIZD!7b>UVM8G5N#%myBeYiZ3tlCoVVNa=V&#K29cyPOw4iog*Ix#X`vp z%2zGQcXB_a;(gX=HoqWC9FU7Q@Ydf@bt9DIcQPpUG51%P{{M>;* z;^_!PvgS^KAOAj}I{k^bLofhbl8plZ*siTkh5mdKq+c`@lIKbuOZ)duMXSaW(L5zI zopE`@2G4q59EBeE@L7wGUGb9Jz^5$+mDS2sO(Uy#?#g258wRgL0W`8w>mxKa^6JU0 z)Yj#bz(Lt=-bPvXSi)VV$yj(OhS6W9%w^$6q_4SSRPp z?{1WlJ#c9xIrbg%@$9qniTq$?4F~eD8}YAP%d{-jpW`{vtW*Q#j(fV^k2WyoT~`&d zXf2oTlJ>OZsr2C`(%x;&q6DDNiW!QjT&9Nn1apilC#A?8OQ$HjP!o4bKARkK-3@Lh zi3fKkNvLOj|KuuHHx6UA4_ZG|$Grpleh1S%*2+9^!+k7DAy?<6|7sI`VX2lZ3vY{@ zg?IEFik!4|{xvVgvu2<9JV6+;I_Pe0>)vI}!mpXtAX%P73#>;&Qsb2T5-JQ12LZ_xs9bfj!z;G>@6|^pGG6BGOe(UqWuXp?huC3nN}Z zP1T>8QM~Q9T0Z$YE>+a0+GYc!W{vIJZgS|J-D%rJk`JG(I-4O00? z>859SPssQWx$#@Cs=U`01ZN$oJ9oDlGz5wCfX zk4%7(n>9x`kj!}bl|X#(>E`LEi!34H@SDBDPYpwli;y@L`N+4EK|RmfZER}qO=lE%HGDASDJsy7 z+@>uD?BHOKNWt09g{@38#m@;bnT-dPjTQR9tuHMHzVZGxnQVl5`GvB5kEV~cG z!sikTr^dtcPQUc9&$gAMTT=Dh41%pP6Ky;ADo9#3qV{>F+4r0b$*8E2 z1j_LaY92RNCO<>xE@)v=fl(VK90fUP4F$s*+Z$kylt)YbtfCH_Ddc(e){^B7rxr&M zrCmg(wBMJOza^s`9M_4-NAUhfZdB-wgsdcA&K1mh+S0==4RBe$hcZcaF~7JoUa9d5OJpK&+#?3yN}C`_hrY0alSUXM7&sQ~7oC zb&F#;gOt+y92{T-%SiK)xT?Da=;&#(n#vT9$&49Boj65L*Xfy$t1Mdokptt9$Z ze=KzkF!T%)iRJLhIY-8gY^t;8`8(E$Q5Y0&@5Kw+!i<=eBeBVx4R_1TUUP3jmMONM z*DYv|WkB^VP6+zW=Ae^}{n_JTqOli?cRZ_gO_5y1=>&mZ4Kc2W<2JLABFaEu1UL;& zh9l-_Qnl6UI)zicUcDI+uXzgfJ+f@*9f~;7%a8_=yVg+&CH+v-I85oweJ4k#Oe#W* zhCWxzwp_n-H1GLH!I_`H)*_u^u@lde zvA3sZ1kHUlqa#KFR|nfvLcTxG+cUOKdbGJc_V&FW8QG-iu^Elp$knHfRD=A%n;Nlu zz(M2Au_lxKecDw9JbIC=#qrbAEAbqxtVu7Elk#`IAt0h(V%cv^p!hloOv395c^Ueai%3kp=xGS0*Yg;R9np#ezq%xkUflTp zCm8c!iOpo(Cd@ckDO5lQ5QfW$?VPv(_wlPQEjMF3T1!Iac=ttXlx`mjfNN8p>A9C5 z9wrbUH5w?AL1PCwXt!Y;>LGzq>Z5~`jqJ&7RQ(qLEZjPUeslZ|oQHEI_}_5JDvDkY zi#gT)wutc8DjsT??L&A$$x#nt@@Axe=%a7ztx!HQHr<^Hz65^SJH>uO%@Y}Z~u8G@JIR3SHxYRlemE*N~zZ2`AgG;6)p$d zOWPKX2?%qZUp~$+J7GYO92Nj$auc#fu9C&_}wU`Q>)OtafWJ6Seu88W4Q&9 z9Q=Ph2+u$scv7F}55sD;gZL6~OT%vm?Yd9FwWr-KWi`(by(B8NuL}KMo}IG>(48Uz zHB___@#^*}$YU70{8<*y4+NgXn^%Ru4>6z8LP`b~yO8Ccr^}rNpt!7j}tvgzW zkWyL;TGF2Hpa12|M~20^@6%92dkxhM3ee8we&Ma?l)^srHuTNm&(Vs2CFua0l|G{? zLAxKuwBo*+baZqbqO^u6Zwnkv*=C8LOqTwF-}71!RAKoo4@avxIfxQgkTM7_Xl8YdU)zM`}--PmMl_Ti{9m*LMez4T)^vTXGKz@4S)H4aH*k zKy5`QWYM@U6$1_`rhcH83=hRODb>%o_q$dHq5fxH_J8e-sHs#VUuK&=A$Ona^I|$H z3jW;`bJ$)(dV>PTcs92w-Vq&8Zg?^}+|_*>&<=rB0nGL9Bxe4_&VF|te>7|E@TlQ{ z=)ldCGIw`JWUbZ7c}vbuMg|Z+^ABOjlwC$dvYE2cDXh>XGQ8so{=}|;-8T->tJD-Q zBaDm9#^tutX|U4sQSTtWcpSC47b3$Zi8ZEVLXla ztp>spyJTl%KL5SP?p%QAu7NnUyK1Vn>F6t>DOPaD9kwx!+u5c;rJiDTj%~|fuF}S& z6{>Yabi#9=PCBh%2C;!6ezKWoo{*WoX$O>VUJedM`=i{cE)N}Yh7`NMK^5=$JDouV zcG8ESmsS(5-n1p&7ygh2UCPX!&i`daRp?D3C2cm!o_b!*8(?FPQv^-)c$ z3neKJQ>>|hjIWiaWAB=CU%F#4GXu+CJg(ds}?6y}((Kp{;G>HRxmN%Nub z)0$Dqpw7vdckkX&30i0LceFVy4wiJ>WCkn(|9WUw^TZ9hjUiMGri#Lgt|ZvmZEglj z8GuQT`m%gDkgW@-E$1tl^a&ym za+a34s5_PPZ!S9kZ>0K))rnQ?`6{#(%Gz?JN*RO_Ec0H%SA!8^@P*PhM>`+XR(vM) zNgjeqdzjW41@VcoYj~a4%znd&(lSnCUPpxe^4t>C??V@5KPD`w`(u;M=Srgw!j4X0q8cLdR#5uc?B>5 zNp-EK`D$h4#rkqBl7b~ry@tASK9}2TGV`kV=crU;nhMoT;a$ER$ZaY?^+F*hKcD!} zvT+K6wejHu`BcbsQBn`+(jo1bMhu{nUY$Xm>y#N zvaDP1i6~roqvvm~Hl0P_m*;lId5!0==+vGMn)Dt{Hua%nzUOtHOkLXpJM0@EVi~M|T>ic5?;`nsMH;)MoNz2G970tES;)UL{;&T3!zKLd k4gMbd|0|9EKcsQOvXCH}IqY~kj0L<@6*Uz~