From 294bde2a474c02bcf07c640a1cf3f26917cc7747 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 1 Dec 2023 15:26:47 +0100 Subject: [PATCH 01/35] monte-carlo improvement --- config.default.yaml | 17 ++- config.tutorial.yaml | 17 ++- scripts/monte_carlo.py | 235 ++++++++++++++++++++++++++++++++++++----- 3 files changed, 238 insertions(+), 31 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index dd05a9fe2..6e11e5cf8 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -343,8 +343,21 @@ costs: monte_carlo: options: add_to_snakefile: false - samples: 7 # number of optimizations - sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported + # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html + # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html + # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html + # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html + # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html + # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html + distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" + # [mean, std] for Normal and LogNormal + # [lower_bound, upper_bound] for Uniform + # [lower_bound, midpoint, upper_bound] for Triangle + # [alpha, beta] for Beta + # [shape, scale] for Gamma + distribution_params: [0,1] + samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number + sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported pypsa_standard: # User can add here flexibly more features for the Monte-Carlo sampling. # Given as "key: value" format diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 6aba3e46c..b082cc7f0 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -339,8 +339,21 @@ costs: monte_carlo: options: add_to_snakefile: false - samples: 7 # number of optimizations - sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported + # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html + # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html + # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html + # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html + # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html + # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html + distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" + # [mean, std] for Normal and LogNormal + # [lower_bound, upper_bound] for Uniform + # [lower_bound, midpoint, upper_bound] for Triangle + # [alpha, beta] for Beta + # [shape, scale] for Gamma + distribution_params: [0,1] + samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number + sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported pypsa_standard: # User can add here flexibly more features for the Monte-Carlo sampling. # Given as "key: value" format diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 5eecb276c..1cbcdf169 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -14,15 +14,30 @@ monte_carlo: options: - add_to_snakefile: - samples: - sampling_strategy: + add_to_snakefile: false + # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html + # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html + # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html + # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html + # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html + # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html + distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" + # [mean, std] for Normal and LogNormal + # [lower_bound, upper_bound] for Uniform + # [lower_bound, midpoint, upper_bound] for Triangle + # [alpha, beta] for Beta + # [shape, scale] for Gamma + distribution_params: [0,1] + samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number + sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported pypsa_standard: - - Examples... + # User can add here flexibly more features for the Monte-Carlo sampling. + # Given as "key: value" format + # Key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! + # Value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object loads_t.p_set: [0.9, 1.1] - generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: [0.9, 1.1] - generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]: [0.9, 1.1] + # generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: [0.9, 1.1] + # generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]: [0.9, 1.1] .. seealso:: Documentation of the configuration file ``config.yaml`` at :ref:`_monte_cf` @@ -81,6 +96,8 @@ def monte_carlo_sampling_pydoe2( N_FEATURES, SAMPLES, + DISTRIBUTION, + DISTRIBUTION_PARAMS, criterion=None, iteration=None, random_state=42, @@ -103,38 +120,65 @@ def monte_carlo_sampling_pydoe2( random_state=random_state, correlation_matrix=correlation_matrix, ) + + lh = rescale_distribution(lh, DISTRIBUTION, DISTRIBUTION_PARAMS) discrepancy = qmc.discrepancy(lh) - logger.info(f"Hypercube discrepancy is: {discrepancy}") + logger.info("Discrepancy is:", discrepancy, + " more details in function documentation.") return lh -def monte_carlo_sampling_chaospy(N_FEATURES, SAMPLES, rule="latin_hypercube", seed=42): +def monte_carlo_sampling_chaospy( + N_FEATURES, + SAMPLES, + DISTRIBUTION, + DISTRIBUTION_PARAMS, + rule="latin_hypercube", + seed=42, +): """ Creates Latin Hypercube Sample (LHS) implementation from chaospy. Documentation on Chaospy: https://github.com/clicumu/pyDOE2 (fixes latin_cube errors) Documentation on Chaospy latin-hyper cube (quasi-Monte Carlo method): https://chaospy.readthedocs.io/en/master/user_guide/fundamentals/quasi_random_samples.html#Quasi-random-samples + """ - # Generate a Nfeatures-dimensional latin hypercube varying between 0 and 1: - N_FEATURES = "chaospy.Uniform(0, 1), " * N_FEATURES - uniform_cube = eval( + import chaospy + from scipy.stats import qmc + from sklearn.preprocessing import MinMaxScaler + + params = tuple(DISTRIBUTION_PARAMS) + # generate a Nfeatures-dimensional latin hypercube varying between 0 and 1: + N_FEATURES = f"chaospy.{DISTRIBUTION}{params}, " * N_FEATURES + cube = eval( f"chaospy.J({N_FEATURES})" ) # writes Nfeatures times the chaospy.uniform... command) - lh = np.atleast_2d(uniform_cube.sample(SAMPLES, rule=rule, seed=seed)).T + lh = cube.sample(SAMPLES, rule=rule, seed=seed).T + + # to check the discrepancy of the samples, the values needs to be from 0-1 + mm = MinMaxScaler(feature_range=(0, 1), clip=True) + lh = mm.fit_transform(lh) + discrepancy = qmc.discrepancy(lh) - logger.info(f"Hypercube discrepancy is: {discrepancy}") + logger.info("Discrepancy is:", discrepancy, + " more details in function documentation.") return lh def monte_carlo_sampling_scipy( - N_FEATURES, SAMPLES, centered=False, strength=2, optimization=None, seed=42 + N_FEATURES, + SAMPLES, + DISTRIBUTION, + DISTRIBUTION_PARAMS, + centered=False, + strength=2, + optimization=None, + seed=42, ): """ - Creates Latin Hypercube Sample (LHS) implementation from SciPy with various - options: - + Creates Latin Hypercube Sample (LHS) implementation from SciPy with various options: - Center the point within the multi-dimensional grid, centered=True - optimization scheme, optimization="random-cd" - strength=1, classical LHS @@ -147,6 +191,8 @@ def monte_carlo_sampling_scipy( Documentation for Latin Hypercube: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.qmc.LatinHypercube.html#scipy.stats.qmc.LatinHypercube Orthogonal LHS is better than basic LHS: https://github.com/scipy/scipy/pull/14546/files, https://en.wikipedia.org/wiki/Latin_hypercube_sampling """ + from scipy.stats import qmc + sampler = qmc.LatinHypercube( d=N_FEATURES, centered=centered, @@ -154,13 +200,135 @@ def monte_carlo_sampling_scipy( optimization=optimization, seed=seed, ) + lh = sampler.random(n=SAMPLES) + + lh = rescale_distribution(lh, DISTRIBUTION, DISTRIBUTION_PARAMS) discrepancy = qmc.discrepancy(lh) - logger.info(f"Hypercube discrepancy is: {discrepancy}") + logger.info("Discrepancy is:", discrepancy, + " more details in function documentation.") return lh +def rescale_distribution(latin_hypercube: np.array, distribution: str, + distribution_params: list): + """ + Rescales a Latin hypercube sampling (LHS) using specified distribution parameters. + + Parameters: + - latin_hypercube (np.array): The Latin hypercube sampling to be rescaled. + - distribution (str): The target distribution for rescaling. Supported options: + "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma". + - distribution_params (list): Parameters specific to the chosen distribution. + For example, for Normal distribution, it should be [mean, std]. + + Returns: + - np.array: Rescaled Latin hypercube sampling. + + Note: + - The function supports rescaling for Uniform, Normal, LogNormal, Triangle, Beta, and Gamma distributions. + - The rescaled samples will have values in the range [0, 1]. + """ + from scipy.stats import norm, lognorm, beta, gamma, triang + from sklearn.preprocessing import MinMaxScaler + + if distribution == "Uniform": + pass + elif distribution == "Normal": + mean, std = distribution_params + latin_hypercube = norm.ppf(latin_hypercube, mean, std) + elif distribution == "LogNormal": + mean, std = distribution_params + latin_hypercube = lognorm.ppf(latin_hypercube, s=0.90) + elif distribution == "Triangle": + tri_mean = np.mean(distribution_params) + latin_hypercube = triang.ppf(latin_hypercube, tri_mean) + elif distribution == "Beta": + a, b = distribution_params + latin_hypercube = beta.ppf(latin_hypercube, a, b) + elif distribution == "Gamma": + shape, scale = distribution_params + latin_hypercube = gamma.ppf(latin_hypercube, shape, scale) + + # samples space needs to be from 0 to 1 + mm = MinMaxScaler(feature_range=(0, 1), clip=True) + latin_hypercube = mm.fit_transform(latin_hypercube) + + return latin_hypercube + + +def validate_parameters( + sampling_strategy: str, samples: int, distribution: str, distribution_params: list +) -> None: + """ + Validates the parameters for a given probability distribution. + Inputs from user through the config file needs to be validated before proceeding to perform monte-carlo simulations. + + Parameters: + - sampling_strategy: str + The chosen sampling strategy from chaospy, scipy and pydoe2 + - samples: int + The number of samples to generate for the simulation + - distribution: str + The name of the probability distribution. + - distribution_params: list + The parameters associated with the probability distribution. + + Raises: + - ValueError: If the parameters are invalid for the specified distribution. + """ + + valid_strategy = ["chaospy", "scipy", "pydoe2"] + valid_distribution = ["Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma"] + + # verifying samples and distribution_params + if samples is None: + raise ValueError(f"assign a value to samples") + elif type(samples) is float or type(samples) is str: + raise ValueError(f"samples must be an integer") + elif distribution_params is None or len(distribution_params) == 0: + raise ValueError(f"assign a list of parameters to distribution_params") + + # verify sampling strategy + if sampling_strategy not in valid_strategy: + raise ValueError( + f" Invalid sampling strategy: {sampling_strategy}. Choose from {valid_strategy}" + ) + + # verify distribution + if distribution not in valid_distribution: + raise ValueError( + f"Unsupported Distribution : {distribution}. Choose from {valid_distribution}" + ) + + # special case handling for Triangle distribution + if distribution == "Triangle": + if len(distribution_params) == 2: + print( + f"{distribution} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" + ) + # use the mean as the middle value + distribution_params.insert(1, np.mean(distribution_params)) + elif len(distribution_params) != 3: + raise ValueError( + f"{distribution} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" + ) + + if distribution in ["Normal", "LogNormal", "Uniform", "Beta", "Gamma"]: + if len(distribution_params) != 2: + raise ValueError(f"{distribution} distribution must have 2 parameters") + + # handling having 0 as values in Beta and Gamma + if distribution in ["Beta", "Gamma"]: + if np.min(distribution_params) <= 0: + raise ValueError( + f"{distribution} distribution cannot have values lower than zero in parameters" + ) + + return None + + if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -180,18 +348,25 @@ def monte_carlo_sampling_scipy( ### SCENARIO INPUTS ### MONTE_CARLO_PYPSA_FEATURES = { - k: v for k, v in monte_carlo_config["pypsa_standard"].items() if v + k: v + for k, v in monte_carlo_config["pypsa_standard"].items() if v } # removes key value pairs with empty value e.g. [] MONTE_CARLO_OPTIONS = monte_carlo_config["options"] L_BOUNDS = [item[0] for item in MONTE_CARLO_PYPSA_FEATURES.values()] U_BOUNDS = [item[1] for item in MONTE_CARLO_PYPSA_FEATURES.values()] - N_FEATURES = len( - MONTE_CARLO_PYPSA_FEATURES - ) # only counts features when specified in config + N_FEATURES = len(MONTE_CARLO_PYPSA_FEATURES + ) # only counts features when specified in config SAMPLES = MONTE_CARLO_OPTIONS.get( "samples" ) # TODO: What is the optimal sampling? Fabian Neumann answered that in "Broad ranges" paper SAMPLING_STRATEGY = MONTE_CARLO_OPTIONS.get("sampling_strategy") + DISTRIBUTION = MONTE_CARLO_OPTIONS.get("distribution") + DISTRIBUTION_PARAMS = MONTE_CARLO_OPTIONS.get("distribution_params") + + ### PARAMETERS VALIDATION + # validates the parameters supplied from config file + validate_parameters(SAMPLING_STRATEGY, SAMPLES, DISTRIBUTION, + DISTRIBUTION_PARAMS) ### SCENARIO CREATION / SAMPLING STRATEGY ### @@ -199,6 +374,8 @@ def monte_carlo_sampling_scipy( lh = monte_carlo_sampling_pydoe2( N_FEATURES, SAMPLES, + DISTRIBUTION, + DISTRIBUTION_PARAMS, criterion=None, iteration=None, random_state=42, @@ -208,6 +385,8 @@ def monte_carlo_sampling_scipy( lh = monte_carlo_sampling_scipy( N_FEATURES, SAMPLES, + DISTRIBUTION, + DISTRIBUTION_PARAMS, centered=False, strength=2, optimization=None, @@ -217,6 +396,8 @@ def monte_carlo_sampling_scipy( lh = monte_carlo_sampling_chaospy( N_FEATURES, SAMPLES, + DISTRIBUTION, + DISTRIBUTION_PARAMS, rule="latin_hypercube", seed=42, ) @@ -235,13 +416,13 @@ def monte_carlo_sampling_scipy( # i, j interaction number to pick values of experimental setup # Example: n.loads_t.p_set = network.loads_t.p_set = .loads_t.p_set * lh_scaled[0,0] exec(f"n.{k} = n.{k} * {lh_scaled[i,j]}") - logger.info(f"Scaled n.{k} by factor {lh_scaled[i,j]} in the {i} scenario") + logger.info( + f"Scaled n.{k} by factor {lh_scaled[i,j]} in the {i} scenario") j = j + 1 ### EXPORT AND METADATA # - latin_hypercube_dict = ( - pd.DataFrame(lh_scaled).rename_axis("Nruns").add_suffix("_feature") - ).to_dict() + latin_hypercube_dict = (pd.DataFrame(lh_scaled).rename_axis( + "Nruns").add_suffix("_feature")).to_dict() n.meta.update(latin_hypercube_dict) n.export_to_netcdf(snakemake.output[0]) From 9a0ddac3574b177d6336a628966d2bed2fbe586a Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 1 Dec 2023 15:31:50 +0100 Subject: [PATCH 02/35] reorganize import --- scripts/monte_carlo.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 1cbcdf169..e96dc1974 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -87,7 +87,9 @@ import pypsa from _helpers import configure_logging, create_logger from pyDOE2 import lhs -from scipy.stats import qmc +import chaospy +from scipy.stats import qmc, norm, lognorm, beta, gamma, triang +from sklearn.preprocessing import MinMaxScaler from solve_network import * logger = create_logger(__name__) @@ -142,11 +144,7 @@ def monte_carlo_sampling_chaospy( Documentation on Chaospy: https://github.com/clicumu/pyDOE2 (fixes latin_cube errors) Documentation on Chaospy latin-hyper cube (quasi-Monte Carlo method): https://chaospy.readthedocs.io/en/master/user_guide/fundamentals/quasi_random_samples.html#Quasi-random-samples - """ - import chaospy - from scipy.stats import qmc - from sklearn.preprocessing import MinMaxScaler params = tuple(DISTRIBUTION_PARAMS) # generate a Nfeatures-dimensional latin hypercube varying between 0 and 1: @@ -191,7 +189,6 @@ def monte_carlo_sampling_scipy( Documentation for Latin Hypercube: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.qmc.LatinHypercube.html#scipy.stats.qmc.LatinHypercube Orthogonal LHS is better than basic LHS: https://github.com/scipy/scipy/pull/14546/files, https://en.wikipedia.org/wiki/Latin_hypercube_sampling """ - from scipy.stats import qmc sampler = qmc.LatinHypercube( d=N_FEATURES, @@ -230,8 +227,6 @@ def rescale_distribution(latin_hypercube: np.array, distribution: str, - The function supports rescaling for Uniform, Normal, LogNormal, Triangle, Beta, and Gamma distributions. - The rescaled samples will have values in the range [0, 1]. """ - from scipy.stats import norm, lognorm, beta, gamma, triang - from sklearn.preprocessing import MinMaxScaler if distribution == "Uniform": pass From 2e34dee15ad8f66611805c13af9cf7049d65a03f Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 1 Dec 2023 17:22:39 +0100 Subject: [PATCH 03/35] precommit --- config.default.yaml | 4 +-- config.tutorial.yaml | 4 +-- scripts/monte_carlo.py | 62 +++++++++++++++++++++++------------------- 3 files changed, 38 insertions(+), 32 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 6e11e5cf8..d40f02444 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -352,10 +352,10 @@ monte_carlo: distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" # [mean, std] for Normal and LogNormal # [lower_bound, upper_bound] for Uniform - # [lower_bound, midpoint, upper_bound] for Triangle + # [lower_bound, midpoint, upper_bound] for Triangle # [alpha, beta] for Beta # [shape, scale] for Gamma - distribution_params: [0,1] + distribution_params: [0, 1] samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported pypsa_standard: diff --git a/config.tutorial.yaml b/config.tutorial.yaml index b082cc7f0..d98f30a14 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -348,10 +348,10 @@ monte_carlo: distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" # [mean, std] for Normal and LogNormal # [lower_bound, upper_bound] for Uniform - # [lower_bound, midpoint, upper_bound] for Triangle + # [lower_bound, midpoint, upper_bound] for Triangle # [alpha, beta] for Beta # [shape, scale] for Gamma - distribution_params: [0,1] + distribution_params: [0, 1] samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported pypsa_standard: diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index e96dc1974..1cde58c19 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -18,13 +18,13 @@ # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html - # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html + # Triangle: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" # [mean, std] for Normal and LogNormal # [lower_bound, upper_bound] for Uniform - # [lower_bound, midpoint, upper_bound] for Triangle + # [lower_bound, midpoint, upper_bound] for Triangle # [alpha, beta] for Beta # [shape, scale] for Gamma distribution_params: [0,1] @@ -87,8 +87,7 @@ import pypsa from _helpers import configure_logging, create_logger from pyDOE2 import lhs -import chaospy -from scipy.stats import qmc, norm, lognorm, beta, gamma, triang +from scipy.stats import beta, gamma, lognorm, norm, qmc, triang from sklearn.preprocessing import MinMaxScaler from solve_network import * @@ -125,8 +124,9 @@ def monte_carlo_sampling_pydoe2( lh = rescale_distribution(lh, DISTRIBUTION, DISTRIBUTION_PARAMS) discrepancy = qmc.discrepancy(lh) - logger.info("Discrepancy is:", discrepancy, - " more details in function documentation.") + logger.info( + "Discrepancy is:", discrepancy, " more details in function documentation." + ) return lh @@ -159,8 +159,9 @@ def monte_carlo_sampling_chaospy( lh = mm.fit_transform(lh) discrepancy = qmc.discrepancy(lh) - logger.info("Discrepancy is:", discrepancy, - " more details in function documentation.") + logger.info( + "Discrepancy is:", discrepancy, " more details in function documentation." + ) return lh @@ -176,7 +177,9 @@ def monte_carlo_sampling_scipy( seed=42, ): """ - Creates Latin Hypercube Sample (LHS) implementation from SciPy with various options: + Creates Latin Hypercube Sample (LHS) implementation from SciPy with various + options: + - Center the point within the multi-dimensional grid, centered=True - optimization scheme, optimization="random-cd" - strength=1, classical LHS @@ -202,22 +205,25 @@ def monte_carlo_sampling_scipy( lh = rescale_distribution(lh, DISTRIBUTION, DISTRIBUTION_PARAMS) discrepancy = qmc.discrepancy(lh) - logger.info("Discrepancy is:", discrepancy, - " more details in function documentation.") + logger.info( + "Discrepancy is:", discrepancy, " more details in function documentation." + ) return lh -def rescale_distribution(latin_hypercube: np.array, distribution: str, - distribution_params: list): +def rescale_distribution( + latin_hypercube: np.array, distribution: str, distribution_params: list +): """ - Rescales a Latin hypercube sampling (LHS) using specified distribution parameters. + Rescales a Latin hypercube sampling (LHS) using specified distribution + parameters. Parameters: - latin_hypercube (np.array): The Latin hypercube sampling to be rescaled. - - distribution (str): The target distribution for rescaling. Supported options: + - distribution (str): The target distribution for rescaling. Supported options: "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma". - - distribution_params (list): Parameters specific to the chosen distribution. + - distribution_params (list): Parameters specific to the chosen distribution. For example, for Normal distribution, it should be [mean, std]. Returns: @@ -257,8 +263,9 @@ def validate_parameters( sampling_strategy: str, samples: int, distribution: str, distribution_params: list ) -> None: """ - Validates the parameters for a given probability distribution. - Inputs from user through the config file needs to be validated before proceeding to perform monte-carlo simulations. + Validates the parameters for a given probability distribution. Inputs from + user through the config file needs to be validated before proceeding to + perform monte-carlo simulations. Parameters: - sampling_strategy: str @@ -343,14 +350,14 @@ def validate_parameters( ### SCENARIO INPUTS ### MONTE_CARLO_PYPSA_FEATURES = { - k: v - for k, v in monte_carlo_config["pypsa_standard"].items() if v + k: v for k, v in monte_carlo_config["pypsa_standard"].items() if v } # removes key value pairs with empty value e.g. [] MONTE_CARLO_OPTIONS = monte_carlo_config["options"] L_BOUNDS = [item[0] for item in MONTE_CARLO_PYPSA_FEATURES.values()] U_BOUNDS = [item[1] for item in MONTE_CARLO_PYPSA_FEATURES.values()] - N_FEATURES = len(MONTE_CARLO_PYPSA_FEATURES - ) # only counts features when specified in config + N_FEATURES = len( + MONTE_CARLO_PYPSA_FEATURES + ) # only counts features when specified in config SAMPLES = MONTE_CARLO_OPTIONS.get( "samples" ) # TODO: What is the optimal sampling? Fabian Neumann answered that in "Broad ranges" paper @@ -360,8 +367,7 @@ def validate_parameters( ### PARAMETERS VALIDATION # validates the parameters supplied from config file - validate_parameters(SAMPLING_STRATEGY, SAMPLES, DISTRIBUTION, - DISTRIBUTION_PARAMS) + validate_parameters(SAMPLING_STRATEGY, SAMPLES, DISTRIBUTION, DISTRIBUTION_PARAMS) ### SCENARIO CREATION / SAMPLING STRATEGY ### @@ -411,13 +417,13 @@ def validate_parameters( # i, j interaction number to pick values of experimental setup # Example: n.loads_t.p_set = network.loads_t.p_set = .loads_t.p_set * lh_scaled[0,0] exec(f"n.{k} = n.{k} * {lh_scaled[i,j]}") - logger.info( - f"Scaled n.{k} by factor {lh_scaled[i,j]} in the {i} scenario") + logger.info(f"Scaled n.{k} by factor {lh_scaled[i,j]} in the {i} scenario") j = j + 1 ### EXPORT AND METADATA # - latin_hypercube_dict = (pd.DataFrame(lh_scaled).rename_axis( - "Nruns").add_suffix("_feature")).to_dict() + latin_hypercube_dict = ( + pd.DataFrame(lh_scaled).rename_axis("Nruns").add_suffix("_feature") + ).to_dict() n.meta.update(latin_hypercube_dict) n.export_to_netcdf(snakemake.output[0]) From a717cf25b548730cb51cc652fba5db047654d3bd Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 12 Dec 2023 20:00:53 +0100 Subject: [PATCH 04/35] revised implementation based on input --- config.default.yaml | 36 +++++---- config.tutorial.yaml | 54 +++++++------ scripts/monte_carlo.py | 171 ++++++++++++++++++++++------------------- 3 files changed, 145 insertions(+), 116 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index d40f02444..56484e313 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -341,31 +341,39 @@ costs: monte_carlo: - options: - add_to_snakefile: false # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" # [mean, std] for Normal and LogNormal # [lower_bound, upper_bound] for Uniform # [lower_bound, midpoint, upper_bound] for Triangle # [alpha, beta] for Beta # [shape, scale] for Gamma - distribution_params: [0, 1] - samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number - sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported - pypsa_standard: - # User can add here flexibly more features for the Monte-Carlo sampling. - # Given as "key: value" format - # Key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! - # Value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object - loads_t.p_set: [0.9, 1.1] - # generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: [0.9, 1.1] - # generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]: [0.9, 1.1] + # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max(): [-3000MW 0MW]"" + # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] + # ... user can add flexibly more features for the Monte-Carlo sampling + # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! + # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object + options: + add_to_snakefile: false + samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number + sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported + uncertainties: + loads_t.p_set: + scale: [0.9, 1.1] + type: Uniform + args: [0, 1] + generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: + scale: [0.9, 1.1] + type: LogNormal + args: [0, 2] + generators_t.p_min_pu.loc[:, n.generators.carrier == "wind"]: + scale: [0.5, 1.2] + type: Beta + args: [0.5, 2] solving: diff --git a/config.tutorial.yaml b/config.tutorial.yaml index d98f30a14..6f522d38c 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -337,31 +337,39 @@ costs: monte_carlo: + # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html + # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html + # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html + # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html + # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html + # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html + # [mean, std] for Normal and LogNormal + # [lower_bound, upper_bound] for Uniform + # [lower_bound, midpoint, upper_bound] for Triangle + # [alpha, beta] for Beta + # [shape, scale] for Gamma + # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max(): [-3000MW 0MW]"" + # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] + # ... user can add flexibly more features for the Monte-Carlo sampling + # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! + # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object options: add_to_snakefile: false - # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html - # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html - # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html - # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html - # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html - # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" - # [mean, std] for Normal and LogNormal - # [lower_bound, upper_bound] for Uniform - # [lower_bound, midpoint, upper_bound] for Triangle - # [alpha, beta] for Beta - # [shape, scale] for Gamma - distribution_params: [0, 1] - samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number - sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported - pypsa_standard: - # User can add here flexibly more features for the Monte-Carlo sampling. - # Given as "key: value" format - # Key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! - # Value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object - loads_t.p_set: [0.9, 1.1] - # generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: [0.9, 1.1] - # generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]: [0.9, 1.1] + samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number + sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported + uncertainties: + loads_t.p_set: + scale: [0.9, 1.1] + type: Uniform + args: [0, 1] + generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: + scale: [0.9, 1.1] + type: LogNormal + args: [0, 2] + generators_t.p_min_pu.loc[:, n.generators.carrier == "wind"]: + scale: [0.5, 1.2] + type: Beta + args: [0.5, 2] solving: diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 1cde58c19..9c5d00108 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -97,8 +97,7 @@ def monte_carlo_sampling_pydoe2( N_FEATURES, SAMPLES, - DISTRIBUTION, - DISTRIBUTION_PARAMS, + uncertainties_values, criterion=None, iteration=None, random_state=42, @@ -111,6 +110,8 @@ def monte_carlo_sampling_pydoe2( Adapted from Disspaset: https://github.com/energy-modelling-toolkit/Dispa-SET/blob/master/scripts/build_and_run_hypercube.py Documentation on PyDOE2: https://github.com/clicumu/pyDOE2 (fixes latin_cube errors) """ + from pyDOE2 import lhs + from scipy.stats import qmc # Generate a Nfeatures-dimensional latin hypercube varying between 0 and 1: lh = lhs( @@ -122,7 +123,7 @@ def monte_carlo_sampling_pydoe2( correlation_matrix=correlation_matrix, ) - lh = rescale_distribution(lh, DISTRIBUTION, DISTRIBUTION_PARAMS) + lh = rescale_distribution(lh, uncertainties_values) discrepancy = qmc.discrepancy(lh) logger.info( "Discrepancy is:", discrepancy, " more details in function documentation." @@ -134,8 +135,7 @@ def monte_carlo_sampling_pydoe2( def monte_carlo_sampling_chaospy( N_FEATURES, SAMPLES, - DISTRIBUTION, - DISTRIBUTION_PARAMS, + uncertainties_values, rule="latin_hypercube", seed=42, ): @@ -145,19 +145,17 @@ def monte_carlo_sampling_chaospy( Documentation on Chaospy: https://github.com/clicumu/pyDOE2 (fixes latin_cube errors) Documentation on Chaospy latin-hyper cube (quasi-Monte Carlo method): https://chaospy.readthedocs.io/en/master/user_guide/fundamentals/quasi_random_samples.html#Quasi-random-samples """ + import chaospy + from scipy.stats import qmc - params = tuple(DISTRIBUTION_PARAMS) # generate a Nfeatures-dimensional latin hypercube varying between 0 and 1: - N_FEATURES = f"chaospy.{DISTRIBUTION}{params}, " * N_FEATURES - cube = eval( + N_FEATURES = "chaospy.Uniform(0, 1), " * N_FEATURES + uniform_cube = eval( f"chaospy.J({N_FEATURES})" ) # writes Nfeatures times the chaospy.uniform... command) - lh = cube.sample(SAMPLES, rule=rule, seed=seed).T - - # to check the discrepancy of the samples, the values needs to be from 0-1 - mm = MinMaxScaler(feature_range=(0, 1), clip=True) - lh = mm.fit_transform(lh) + lh = uniform_cube.sample(SAMPLES, rule=rule, seed=seed).T + lh = rescale_distribution(lh, uncertainties_values) discrepancy = qmc.discrepancy(lh) logger.info( "Discrepancy is:", discrepancy, " more details in function documentation." @@ -169,8 +167,7 @@ def monte_carlo_sampling_chaospy( def monte_carlo_sampling_scipy( N_FEATURES, SAMPLES, - DISTRIBUTION, - DISTRIBUTION_PARAMS, + uncertainties_values, centered=False, strength=2, optimization=None, @@ -192,6 +189,7 @@ def monte_carlo_sampling_scipy( Documentation for Latin Hypercube: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.qmc.LatinHypercube.html#scipy.stats.qmc.LatinHypercube Orthogonal LHS is better than basic LHS: https://github.com/scipy/scipy/pull/14546/files, https://en.wikipedia.org/wiki/Latin_hypercube_sampling """ + from scipy.stats import qmc sampler = qmc.LatinHypercube( d=N_FEATURES, @@ -203,7 +201,7 @@ def monte_carlo_sampling_scipy( lh = sampler.random(n=SAMPLES) - lh = rescale_distribution(lh, DISTRIBUTION, DISTRIBUTION_PARAMS) + lh = rescale_distribution(lh, uncertainties_values) discrepancy = qmc.discrepancy(lh) logger.info( "Discrepancy is:", discrepancy, " more details in function documentation." @@ -213,44 +211,57 @@ def monte_carlo_sampling_scipy( def rescale_distribution( - latin_hypercube: np.array, distribution: str, distribution_params: list -): + latin_hypercube: np.array, uncertainties_values: dict +) -> np.array: """ Rescales a Latin hypercube sampling (LHS) using specified distribution parameters. Parameters: - latin_hypercube (np.array): The Latin hypercube sampling to be rescaled. - - distribution (str): The target distribution for rescaling. Supported options: - "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma". - - distribution_params (list): Parameters specific to the chosen distribution. - For example, for Normal distribution, it should be [mean, std]. + - uncertainties_values (list): List of dictionaries containing distribution information. + Each dictionary should have 'type' key specifying the distribution type and 'args' key + containing parameters specific to the chosen distribution. Returns: - - np.array: Rescaled Latin hypercube sampling. + - np.array: Rescaled Latin hypercube sampling with values in the range [0, 1]. + + Supported Distributions: + - "Uniform": No rescaling applied. + - "Normal": Rescaled using the inverse of the normal distribution function with specified mean and std. + - "LogNormal": Rescaled using the inverse of the log-normal distribution function with specified mean and std. + - "Triangle": Rescaled using the inverse of the triangular distribution function with mean calculated from given parameters. + - "Beta": Rescaled using the inverse of the beta distribution function with specified shape parameters. + - "Gamma": Rescaled using the inverse of the gamma distribution function with specified shape and scale parameters. Note: - The function supports rescaling for Uniform, Normal, LogNormal, Triangle, Beta, and Gamma distributions. - The rescaled samples will have values in the range [0, 1]. """ - - if distribution == "Uniform": - pass - elif distribution == "Normal": - mean, std = distribution_params - latin_hypercube = norm.ppf(latin_hypercube, mean, std) - elif distribution == "LogNormal": - mean, std = distribution_params - latin_hypercube = lognorm.ppf(latin_hypercube, s=0.90) - elif distribution == "Triangle": - tri_mean = np.mean(distribution_params) - latin_hypercube = triang.ppf(latin_hypercube, tri_mean) - elif distribution == "Beta": - a, b = distribution_params - latin_hypercube = beta.ppf(latin_hypercube, a, b) - elif distribution == "Gamma": - shape, scale = distribution_params - latin_hypercube = gamma.ppf(latin_hypercube, shape, scale) + from scipy.stats import beta, gamma, lognorm, norm, triang + from sklearn.preprocessing import MinMaxScaler + + for idx, value in enumerate(uncertainties_values): + dist = value.get("type") + params = value.get("args") + + if dist == "Uniform": + pass + elif dist == "Normal": + mean, std = params + latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) + elif dist == "LogNormal": + mean, std = params + latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=0.90) + elif dist == "Triangle": + tri_mean = np.mean(params) + latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) + elif dist == "Beta": + a, b = params + latin_hypercube[:, idx] = beta.ppf(latin_hypercube[:, idx], a, b) + elif dist == "Gamma": + shape, scale = params + latin_hypercube[:, idx] = gamma.ppf(latin_hypercube[:, idx], shape, scale) # samples space needs to be from 0 to 1 mm = MinMaxScaler(feature_range=(0, 1), clip=True) @@ -260,7 +271,7 @@ def rescale_distribution( def validate_parameters( - sampling_strategy: str, samples: int, distribution: str, distribution_params: list + sampling_strategy: str, samples: int, uncertainties_values: dict ) -> None: """ Validates the parameters for a given probability distribution. Inputs from @@ -289,8 +300,8 @@ def validate_parameters( raise ValueError(f"assign a value to samples") elif type(samples) is float or type(samples) is str: raise ValueError(f"samples must be an integer") - elif distribution_params is None or len(distribution_params) == 0: - raise ValueError(f"assign a list of parameters to distribution_params") + # elif distribution_params is None or len(distribution_params) == 0: + # raise ValueError(f"assign a list of parameters to distribution_params") # verify sampling strategy if sampling_strategy not in valid_strategy: @@ -298,36 +309,42 @@ def validate_parameters( f" Invalid sampling strategy: {sampling_strategy}. Choose from {valid_strategy}" ) - # verify distribution - if distribution not in valid_distribution: - raise ValueError( - f"Unsupported Distribution : {distribution}. Choose from {valid_distribution}" - ) - - # special case handling for Triangle distribution - if distribution == "Triangle": - if len(distribution_params) == 2: - print( - f"{distribution} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" - ) - # use the mean as the middle value - distribution_params.insert(1, np.mean(distribution_params)) - elif len(distribution_params) != 3: - raise ValueError( - f"{distribution} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" - ) + for idx, value in enumerate(uncertainties_values): + dist_type = value.get("type") + param = value.get("args") + scale = value.get("scale") - if distribution in ["Normal", "LogNormal", "Uniform", "Beta", "Gamma"]: - if len(distribution_params) != 2: - raise ValueError(f"{distribution} distribution must have 2 parameters") + if dist_type is None or len(param) == 0 or len(scale) == 0: + raise ValueError(f"assign a list of parameters to distribution_params") - # handling having 0 as values in Beta and Gamma - if distribution in ["Beta", "Gamma"]: - if np.min(distribution_params) <= 0: + if dist_type not in valid_distribution: raise ValueError( - f"{distribution} distribution cannot have values lower than zero in parameters" + f"Unsupported Distribution : {dist_type}. Choose from {valid_distribution}" ) + if dist_type == "Triangle": + if len(param) == 2: + print( + f"{dist_type} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" + ) + # use the mean as the middle value + param.insert(1, np.mean(param)) + elif len(param) != 3: + raise ValueError( + f"{dist_type} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" + ) + + if dist_type in ["Normal", "LogNormal", "Uniform", "Beta", "Gamma"]: + if len(param) != 2: + raise ValueError(f"{dist_type} distribution must have 2 parameters") + + # handling having 0 as values in Beta and Gamma + if dist_type in ["Beta", "Gamma"]: + if np.min(param) <= 0: + raise ValueError( + f"{dist_type} distribution cannot have values lower than zero in parameters" + ) + return None @@ -350,7 +367,7 @@ def validate_parameters( ### SCENARIO INPUTS ### MONTE_CARLO_PYPSA_FEATURES = { - k: v for k, v in monte_carlo_config["pypsa_standard"].items() if v + k: v.get("scale") for k, v in monte_carlo_config["uncertainties"].items() if v } # removes key value pairs with empty value e.g. [] MONTE_CARLO_OPTIONS = monte_carlo_config["options"] L_BOUNDS = [item[0] for item in MONTE_CARLO_PYPSA_FEATURES.values()] @@ -362,12 +379,11 @@ def validate_parameters( "samples" ) # TODO: What is the optimal sampling? Fabian Neumann answered that in "Broad ranges" paper SAMPLING_STRATEGY = MONTE_CARLO_OPTIONS.get("sampling_strategy") - DISTRIBUTION = MONTE_CARLO_OPTIONS.get("distribution") - DISTRIBUTION_PARAMS = MONTE_CARLO_OPTIONS.get("distribution_params") + uncertainties_values = monte_carlo_config["uncertainties"].values() ### PARAMETERS VALIDATION # validates the parameters supplied from config file - validate_parameters(SAMPLING_STRATEGY, SAMPLES, DISTRIBUTION, DISTRIBUTION_PARAMS) + validate_parameters(SAMPLING_STRATEGY, SAMPLES, uncertainties_values) ### SCENARIO CREATION / SAMPLING STRATEGY ### @@ -375,8 +391,7 @@ def validate_parameters( lh = monte_carlo_sampling_pydoe2( N_FEATURES, SAMPLES, - DISTRIBUTION, - DISTRIBUTION_PARAMS, + uncertainties_values, criterion=None, iteration=None, random_state=42, @@ -386,8 +401,7 @@ def validate_parameters( lh = monte_carlo_sampling_scipy( N_FEATURES, SAMPLES, - DISTRIBUTION, - DISTRIBUTION_PARAMS, + uncertainties_values, centered=False, strength=2, optimization=None, @@ -397,8 +411,7 @@ def validate_parameters( lh = monte_carlo_sampling_chaospy( N_FEATURES, SAMPLES, - DISTRIBUTION, - DISTRIBUTION_PARAMS, + uncertainties_values, rule="latin_hypercube", seed=42, ) From feca2a4f8a1f8fce524d2cb5ecbfccf696317cf5 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 12 Dec 2023 20:07:18 +0100 Subject: [PATCH 05/35] removed unused comments --- scripts/monte_carlo.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 9c5d00108..fd637b325 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -300,8 +300,6 @@ def validate_parameters( raise ValueError(f"assign a value to samples") elif type(samples) is float or type(samples) is str: raise ValueError(f"samples must be an integer") - # elif distribution_params is None or len(distribution_params) == 0: - # raise ValueError(f"assign a list of parameters to distribution_params") # verify sampling strategy if sampling_strategy not in valid_strategy: From 3585bd9c79d08b3b1b3d1e55b6e1598fa8527b21 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 15 Dec 2023 22:53:33 +0100 Subject: [PATCH 06/35] reviewed implementation on monte-carlo --- config.default.yaml | 7 ++---- config.tutorial.yaml | 3 --- scripts/monte_carlo.py | 54 +++++++++++++++++++----------------------- 3 files changed, 27 insertions(+), 37 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 56484e313..fa880de18 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -359,19 +359,16 @@ monte_carlo: # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object options: add_to_snakefile: false - samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number - sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported + samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number + sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported uncertainties: loads_t.p_set: - scale: [0.9, 1.1] type: Uniform args: [0, 1] generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: - scale: [0.9, 1.1] type: LogNormal args: [0, 2] generators_t.p_min_pu.loc[:, n.generators.carrier == "wind"]: - scale: [0.5, 1.2] type: Beta args: [0.5, 2] diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 6f522d38c..6833a99c0 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -359,15 +359,12 @@ monte_carlo: sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported uncertainties: loads_t.p_set: - scale: [0.9, 1.1] type: Uniform args: [0, 1] generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: - scale: [0.9, 1.1] type: LogNormal args: [0, 2] generators_t.p_min_pu.loc[:, n.generators.carrier == "wind"]: - scale: [0.5, 1.2] type: Beta args: [0.5, 2] diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index fd637b325..5f38aa515 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -210,9 +210,7 @@ def monte_carlo_sampling_scipy( return lh -def rescale_distribution( - latin_hypercube: np.array, uncertainties_values: dict -) -> np.array: +def rescale_distribution(latin_hypercube, uncertainties_values): """ Rescales a Latin hypercube sampling (LHS) using specified distribution parameters. @@ -245,23 +243,26 @@ def rescale_distribution( dist = value.get("type") params = value.get("args") - if dist == "Uniform": - pass - elif dist == "Normal": - mean, std = params - latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) - elif dist == "LogNormal": - mean, std = params - latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=0.90) - elif dist == "Triangle": - tri_mean = np.mean(params) - latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) - elif dist == "Beta": - a, b = params - latin_hypercube[:, idx] = beta.ppf(latin_hypercube[:, idx], a, b) - elif dist == "Gamma": - shape, scale = params - latin_hypercube[:, idx] = gamma.ppf(latin_hypercube[:, idx], shape, scale) + match dist: + case "Uniform": + pass + case "Normal": + mean, std = params + latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) + case "LogNormal": + mean, std = params + latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=0.90) + case "Triangle": + tri_mean = np.mean(params) + latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) + case "Beta": + a, b = params + latin_hypercube[:, idx] = beta.ppf(latin_hypercube[:, idx], a, b) + case "Gamma": + shape, scale = params + latin_hypercube[:, idx] = gamma.ppf( + latin_hypercube[:, idx], shape, scale + ) # samples space needs to be from 0 to 1 mm = MinMaxScaler(feature_range=(0, 1), clip=True) @@ -270,9 +271,7 @@ def rescale_distribution( return latin_hypercube -def validate_parameters( - sampling_strategy: str, samples: int, uncertainties_values: dict -) -> None: +def validate_parameters(sampling_strategy, samples, uncertainties_values): """ Validates the parameters for a given probability distribution. Inputs from user through the config file needs to be validated before proceeding to @@ -368,8 +367,6 @@ def validate_parameters( k: v.get("scale") for k, v in monte_carlo_config["uncertainties"].items() if v } # removes key value pairs with empty value e.g. [] MONTE_CARLO_OPTIONS = monte_carlo_config["options"] - L_BOUNDS = [item[0] for item in MONTE_CARLO_PYPSA_FEATURES.values()] - U_BOUNDS = [item[1] for item in MONTE_CARLO_PYPSA_FEATURES.values()] N_FEATURES = len( MONTE_CARLO_PYPSA_FEATURES ) # only counts features when specified in config @@ -413,7 +410,6 @@ def validate_parameters( rule="latin_hypercube", seed=42, ) - lh_scaled = qmc.scale(lh, L_BOUNDS, U_BOUNDS) ### MONTE-CARLO MODIFICATIONS ### @@ -427,14 +423,14 @@ def validate_parameters( # v is the lower and upper bound [0.8,1.3], that was used for lh_scaled # i, j interaction number to pick values of experimental setup # Example: n.loads_t.p_set = network.loads_t.p_set = .loads_t.p_set * lh_scaled[0,0] - exec(f"n.{k} = n.{k} * {lh_scaled[i,j]}") - logger.info(f"Scaled n.{k} by factor {lh_scaled[i,j]} in the {i} scenario") + exec(f"n.{k} = n.{k} * {lh[i,j]}") + logger.info(f"Scaled n.{k} by factor {lh[i,j]} in the {i} scenario") j = j + 1 ### EXPORT AND METADATA # latin_hypercube_dict = ( - pd.DataFrame(lh_scaled).rename_axis("Nruns").add_suffix("_feature") + pd.DataFrame(lh).rename_axis("Nruns").add_suffix("_feature") ).to_dict() n.meta.update(latin_hypercube_dict) n.export_to_netcdf(snakemake.output[0]) From 24c3de76cfc3bb3a31ea95dc0231ffbff49f8535 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Mon, 18 Dec 2023 01:13:47 +0100 Subject: [PATCH 07/35] removed var scale --- scripts/monte_carlo.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 5f38aa515..0a8dd3079 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -309,9 +309,8 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): for idx, value in enumerate(uncertainties_values): dist_type = value.get("type") param = value.get("args") - scale = value.get("scale") - if dist_type is None or len(param) == 0 or len(scale) == 0: + if dist_type is None or len(param) == 0: raise ValueError(f"assign a list of parameters to distribution_params") if dist_type not in valid_distribution: From e79566a271a9d3fd11d3035b33dc5e987129a827 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Mon, 18 Dec 2023 15:15:50 +0100 Subject: [PATCH 08/35] added seedling option --- config.default.yaml | 1 + config.tutorial.yaml | 1 + scripts/monte_carlo.py | 25 +++++++++++++------------ 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index fa880de18..1d3746d0b 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -361,6 +361,7 @@ monte_carlo: add_to_snakefile: false samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported + seed: 42 # set seedling for reproducibilty uncertainties: loads_t.p_set: type: Uniform diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 6833a99c0..464c61502 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -357,6 +357,7 @@ monte_carlo: add_to_snakefile: false samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported + seed: 42 # set seedling for reproducibilty uncertainties: loads_t.p_set: type: Uniform diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 0a8dd3079..6e980e78f 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -98,10 +98,10 @@ def monte_carlo_sampling_pydoe2( N_FEATURES, SAMPLES, uncertainties_values, + random_state, criterion=None, iteration=None, - random_state=42, - correlation_matrix=None, + correlation_matrix=None ): """ Creates Latin Hypercube Sample (LHS) implementation from PyDOE2 with @@ -136,8 +136,8 @@ def monte_carlo_sampling_chaospy( N_FEATURES, SAMPLES, uncertainties_values, - rule="latin_hypercube", - seed=42, + seed, + rule="latin_hypercube" ): """ Creates Latin Hypercube Sample (LHS) implementation from chaospy. @@ -168,10 +168,10 @@ def monte_carlo_sampling_scipy( N_FEATURES, SAMPLES, uncertainties_values, + seed, centered=False, strength=2, - optimization=None, - seed=42, + optimization=None ): """ Creates Latin Hypercube Sample (LHS) implementation from SciPy with various @@ -374,6 +374,7 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): ) # TODO: What is the optimal sampling? Fabian Neumann answered that in "Broad ranges" paper SAMPLING_STRATEGY = MONTE_CARLO_OPTIONS.get("sampling_strategy") uncertainties_values = monte_carlo_config["uncertainties"].values() + seed = MONTE_CARLO_OPTIONS.get("seed") ### PARAMETERS VALIDATION # validates the parameters supplied from config file @@ -386,28 +387,28 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): N_FEATURES, SAMPLES, uncertainties_values, + random_state=seed, criterion=None, iteration=None, - random_state=42, - correlation_matrix=None, + correlation_matrix=None ) if SAMPLING_STRATEGY == "scipy": lh = monte_carlo_sampling_scipy( N_FEATURES, SAMPLES, uncertainties_values, + seed=seed, centered=False, strength=2, - optimization=None, - seed=42, + optimization=None ) if SAMPLING_STRATEGY == "chaospy": lh = monte_carlo_sampling_chaospy( N_FEATURES, SAMPLES, uncertainties_values, - rule="latin_hypercube", - seed=42, + seed=seed, + rule="latin_hypercube" ) ### MONTE-CARLO MODIFICATIONS From f5b62c74e4048752c3dca554e89f2dc9ce158876 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 19 Dec 2023 00:31:37 +0100 Subject: [PATCH 09/35] precommit --- scripts/monte_carlo.py | 62 ++++++++++++++++++++---------------------- 1 file changed, 29 insertions(+), 33 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 6e980e78f..1af3c8de1 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -101,7 +101,7 @@ def monte_carlo_sampling_pydoe2( random_state, criterion=None, iteration=None, - correlation_matrix=None + correlation_matrix=None, ): """ Creates Latin Hypercube Sample (LHS) implementation from PyDOE2 with @@ -133,11 +133,7 @@ def monte_carlo_sampling_pydoe2( def monte_carlo_sampling_chaospy( - N_FEATURES, - SAMPLES, - uncertainties_values, - seed, - rule="latin_hypercube" + N_FEATURES, SAMPLES, uncertainties_values, seed, rule="latin_hypercube" ): """ Creates Latin Hypercube Sample (LHS) implementation from chaospy. @@ -171,7 +167,7 @@ def monte_carlo_sampling_scipy( seed, centered=False, strength=2, - optimization=None + optimization=None, ): """ Creates Latin Hypercube Sample (LHS) implementation from SciPy with various @@ -236,8 +232,8 @@ def rescale_distribution(latin_hypercube, uncertainties_values): - The function supports rescaling for Uniform, Normal, LogNormal, Triangle, Beta, and Gamma distributions. - The rescaled samples will have values in the range [0, 1]. """ - from scipy.stats import beta, gamma, lognorm, norm, triang - from sklearn.preprocessing import MinMaxScaler + from scipy.stats import beta, gamma, lognorm, norm, qmc, triang + from sklearn.preprocessing import MinMaxScaler, minmax_scale for idx, value in enumerate(uncertainties_values): dist = value.get("type") @@ -245,12 +241,15 @@ def rescale_distribution(latin_hypercube, uncertainties_values): match dist: case "Uniform": - pass + l_bounds, u_bounds = params + latin_hypercube[:, idx] = minmax_scale( + latin_hypercube[:, idx], feature_range=(l_bounds, u_bounds) + ) case "Normal": mean, std = params latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) case "LogNormal": - mean, std = params + shape = params[0] latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=0.90) case "Triangle": tri_mean = np.mean(params) @@ -330,9 +329,10 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): f"{dist_type} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" ) - if dist_type in ["Normal", "LogNormal", "Uniform", "Beta", "Gamma"]: - if len(param) != 2: - raise ValueError(f"{dist_type} distribution must have 2 parameters") + if dist_type in ["Normal", "Uniform", "Beta", "Gamma"] and len(param) != 2: + raise ValueError(f"{dist_type} distribution must have 2 parameters") + elif dist_type == "LogNormal" and len(dist_type) == 1: + raise ValueError(f"{dist_type} must have a single parameter") # handling having 0 as values in Beta and Gamma if dist_type in ["Beta", "Gamma"]: @@ -360,7 +360,7 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): configure_logging(snakemake) monte_carlo_config = snakemake.params.monte_carlo - ### SCENARIO INPUTS + # SCENARIO INPUTS ### MONTE_CARLO_PYPSA_FEATURES = { k: v.get("scale") for k, v in monte_carlo_config["uncertainties"].items() if v @@ -373,45 +373,41 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): "samples" ) # TODO: What is the optimal sampling? Fabian Neumann answered that in "Broad ranges" paper SAMPLING_STRATEGY = MONTE_CARLO_OPTIONS.get("sampling_strategy") - uncertainties_values = monte_carlo_config["uncertainties"].values() - seed = MONTE_CARLO_OPTIONS.get("seed") + UNCERTAINTIES_VALUES = monte_carlo_config["uncertainties"].values() + SEED = MONTE_CARLO_OPTIONS.get("seed") - ### PARAMETERS VALIDATION + # PARAMETERS VALIDATION # validates the parameters supplied from config file - validate_parameters(SAMPLING_STRATEGY, SAMPLES, uncertainties_values) + validate_parameters(SAMPLING_STRATEGY, SAMPLES, UNCERTAINTIES_VALUES) - ### SCENARIO CREATION / SAMPLING STRATEGY + # SCENARIO CREATION / SAMPLING STRATEGY ### if SAMPLING_STRATEGY == "pydoe2": lh = monte_carlo_sampling_pydoe2( N_FEATURES, SAMPLES, - uncertainties_values, - random_state=seed, + UNCERTAINTIES_VALUES, + random_state=SEED, criterion=None, iteration=None, - correlation_matrix=None + correlation_matrix=None, ) if SAMPLING_STRATEGY == "scipy": lh = monte_carlo_sampling_scipy( N_FEATURES, SAMPLES, - uncertainties_values, - seed=seed, + UNCERTAINTIES_VALUES, + seed=SEED, centered=False, strength=2, - optimization=None + optimization=None, ) if SAMPLING_STRATEGY == "chaospy": lh = monte_carlo_sampling_chaospy( - N_FEATURES, - SAMPLES, - uncertainties_values, - seed=seed, - rule="latin_hypercube" + N_FEATURES, SAMPLES, UNCERTAINTIES_VALUES, seed=SEED, rule="latin_hypercube" ) - ### MONTE-CARLO MODIFICATIONS + # MONTE-CARLO MODIFICATIONS ### n = pypsa.Network(snakemake.input[0]) unc_wildcards = snakemake.wildcards[-1] @@ -427,7 +423,7 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): logger.info(f"Scaled n.{k} by factor {lh[i,j]} in the {i} scenario") j = j + 1 - ### EXPORT AND METADATA + # EXPORT AND METADATA # latin_hypercube_dict = ( pd.DataFrame(lh).rename_axis("Nruns").add_suffix("_feature") From a8b7809847668e6f88b7529c5122a4fb84a11726 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 19 Dec 2023 12:16:01 +0100 Subject: [PATCH 10/35] fixed legacy scale --- scripts/monte_carlo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 1af3c8de1..99771beb2 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -250,7 +250,7 @@ def rescale_distribution(latin_hypercube, uncertainties_values): latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) case "LogNormal": shape = params[0] - latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=0.90) + latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=shape) case "Triangle": tri_mean = np.mean(params) latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) @@ -363,7 +363,7 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): # SCENARIO INPUTS ### MONTE_CARLO_PYPSA_FEATURES = { - k: v.get("scale") for k, v in monte_carlo_config["uncertainties"].items() if v + k for k in monte_carlo_config["uncertainties"].keys() if k } # removes key value pairs with empty value e.g. [] MONTE_CARLO_OPTIONS = monte_carlo_config["options"] N_FEATURES = len( @@ -413,7 +413,7 @@ def validate_parameters(sampling_strategy, samples, uncertainties_values): unc_wildcards = snakemake.wildcards[-1] i = int(unc_wildcards[1:]) j = 0 - for k, v in MONTE_CARLO_PYPSA_FEATURES.items(): + for k in MONTE_CARLO_PYPSA_FEATURES: # this loop sets in one scenario each "i" feature assumption # k is the config input key "loads_t.p_set" # v is the lower and upper bound [0.8,1.3], that was used for lh_scaled From 9ee395c077deb9e4253cc99147ef4424a1abea4b Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 19 Dec 2023 13:42:10 +0100 Subject: [PATCH 11/35] added type hinting --- scripts/monte_carlo.py | 48 ++++++++++++++++++++++++------------------ 1 file changed, 28 insertions(+), 20 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 99771beb2..763aa1214 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -95,14 +95,14 @@ def monte_carlo_sampling_pydoe2( - N_FEATURES, - SAMPLES, - uncertainties_values, - random_state, - criterion=None, - iteration=None, - correlation_matrix=None, -): + N_FEATURES: int, + SAMPLES: int, + uncertainties_values: dict, + random_state: int, + criterion: str = None, + iteration: int = None, + correlation_matrix: np.ndarray = None, +) -> np.ndarray: """ Creates Latin Hypercube Sample (LHS) implementation from PyDOE2 with various options. Additionally all "corners" are simulated. @@ -133,8 +133,12 @@ def monte_carlo_sampling_pydoe2( def monte_carlo_sampling_chaospy( - N_FEATURES, SAMPLES, uncertainties_values, seed, rule="latin_hypercube" -): + N_FEATURES: int, + SAMPLES: int, + uncertainties_values: dict, + seed: int, + rule: str = "latin_hypercube", +) -> np.ndarray: """ Creates Latin Hypercube Sample (LHS) implementation from chaospy. @@ -161,14 +165,14 @@ def monte_carlo_sampling_chaospy( def monte_carlo_sampling_scipy( - N_FEATURES, - SAMPLES, - uncertainties_values, - seed, - centered=False, - strength=2, - optimization=None, -): + N_FEATURES: int, + SAMPLES: int, + uncertainties_values: dict, + seed: int, + centered: bool = False, + strength: int = 2, + optimization: str = None, +) -> np.ndarray: """ Creates Latin Hypercube Sample (LHS) implementation from SciPy with various options: @@ -206,7 +210,9 @@ def monte_carlo_sampling_scipy( return lh -def rescale_distribution(latin_hypercube, uncertainties_values): +def rescale_distribution( + latin_hypercube: np.ndarray, uncertainties_values: dict +) -> np.ndarray: """ Rescales a Latin hypercube sampling (LHS) using specified distribution parameters. @@ -270,7 +276,9 @@ def rescale_distribution(latin_hypercube, uncertainties_values): return latin_hypercube -def validate_parameters(sampling_strategy, samples, uncertainties_values): +def validate_parameters( + sampling_strategy: str, samples: int, uncertainties_values: dict +) -> None: """ Validates the parameters for a given probability distribution. Inputs from user through the config file needs to be validated before proceeding to From 37c0b7aab0182dbf232ce7f29bec233448b1a413 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Thu, 21 Dec 2023 00:09:22 +0100 Subject: [PATCH 12/35] added feature to plot distribution --- scripts/monte_carlo.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 763aa1214..834bf1d29 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -85,6 +85,7 @@ import numpy as np import pandas as pd import pypsa +import seaborn as sns from _helpers import configure_logging, create_logger from pyDOE2 import lhs from scipy.stats import beta, gamma, lognorm, norm, qmc, triang @@ -92,6 +93,7 @@ from solve_network import * logger = create_logger(__name__) +sns.set() def monte_carlo_sampling_pydoe2( @@ -370,9 +372,9 @@ def validate_parameters( # SCENARIO INPUTS ### - MONTE_CARLO_PYPSA_FEATURES = { + MONTE_CARLO_PYPSA_FEATURES = [ k for k in monte_carlo_config["uncertainties"].keys() if k - } # removes key value pairs with empty value e.g. [] + ] # removes key value pairs with empty value e.g. [] MONTE_CARLO_OPTIONS = monte_carlo_config["options"] N_FEATURES = len( MONTE_CARLO_PYPSA_FEATURES @@ -415,6 +417,12 @@ def validate_parameters( N_FEATURES, SAMPLES, UNCERTAINTIES_VALUES, seed=SEED, rule="latin_hypercube" ) + # create plot for the rescaled distributions + for idx in range(N_FEATURES): + sns.displot(lh[:, idx]).set( + title=f"{MONTE_CARLO_PYPSA_FEATURES[idx]}" + ).figure.savefig(f"{MONTE_CARLO_PYPSA_FEATURES[idx]}.png", bbox_inches="tight") + # MONTE-CARLO MODIFICATIONS ### n = pypsa.Network(snakemake.input[0]) From 1ef294a74a21da7d5de15b8d74b776d7717e2a22 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Thu, 21 Dec 2023 00:12:12 +0100 Subject: [PATCH 13/35] added new feature to clean rule with monte-carlo --- Snakefile | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/Snakefile b/Snakefile index ce3bd9f13..a82d4ce4a 100644 --- a/Snakefile +++ b/Snakefile @@ -65,9 +65,8 @@ if config["custom_rules"] is not []: include: rule -rule clean: +rule clean_monte: run: - shell("snakemake -j 1 solve_all_networks --delete-all-output") try: shell("snakemake -j 1 solve_all_networks_monte --delete-all-output") except: @@ -75,6 +74,15 @@ rule clean: shell("snakemake -j 1 run_all_scenarios --delete-all-output") +rule clean: + run: + try: + shell("snakemake -j 1 solve_all_networks --delete-all-output") + except: + pass + shell("snakemake -j 1 run_all_scenarios --delete-all-output") + + rule run_tests: output: touch("tests.done"), From 72108afc5ff13a02ed3d9c0b09ce8debe0f48b7a Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 22 Dec 2023 09:13:20 +0100 Subject: [PATCH 14/35] fixed lognormal shape --- config.default.yaml | 9 +++++---- config.tutorial.yaml | 5 +++-- scripts/monte_carlo.py | 23 +++++++++++++++-------- 3 files changed, 23 insertions(+), 14 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 09498c57b..46b0d1d1e 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -25,7 +25,7 @@ enable: custom_rules: [] # Default empty [] or link to custom rule file e.g. ["my_folder/my_rules.smk"] that add rules to Snakefile run: - name: "" # use this to keep track of runs with different settings + name: "NG_BJ" # use this to keep track of runs with different settings shared_cutouts: true # set to true to share the default cutout(s) across runs # Note: value false requires build_cutout to be enabled @@ -347,7 +347,8 @@ monte_carlo: # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - # [mean, std] for Normal and LogNormal + # [mean, std] for Normal + # [shape] for LogNormal # [lower_bound, upper_bound] for Uniform # [lower_bound, midpoint, upper_bound] for Triangle # [alpha, beta] for Beta @@ -358,7 +359,7 @@ monte_carlo: # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object options: - add_to_snakefile: false + add_to_snakefile: true samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty @@ -368,7 +369,7 @@ monte_carlo: args: [0, 1] generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: type: LogNormal - args: [0, 2] + args: [1.5] generators_t.p_min_pu.loc[:, n.generators.carrier == "wind"]: type: Beta args: [0.5, 2] diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 464c61502..d25352274 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -343,7 +343,8 @@ monte_carlo: # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - # [mean, std] for Normal and LogNormal + # [mean, std] for Normal + # [shape] for LogNormal # [lower_bound, upper_bound] for Uniform # [lower_bound, midpoint, upper_bound] for Triangle # [alpha, beta] for Beta @@ -364,7 +365,7 @@ monte_carlo: args: [0, 1] generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: type: LogNormal - args: [0, 2] + args: [1.5] generators_t.p_min_pu.loc[:, n.generators.carrier == "wind"]: type: Beta args: [0.5, 2] diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 834bf1d29..6ddd21b09 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -255,16 +255,20 @@ def rescale_distribution( ) case "Normal": mean, std = params - latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) + latin_hypercube[:, idx] = norm.ppf( + latin_hypercube[:, idx], mean, std) case "LogNormal": shape = params[0] - latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=shape) + latin_hypercube[:, idx] = lognorm.ppf( + latin_hypercube[:, idx], s=shape) case "Triangle": tri_mean = np.mean(params) - latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) + latin_hypercube[:, idx] = triang.ppf( + latin_hypercube[:, idx], tri_mean) case "Beta": a, b = params - latin_hypercube[:, idx] = beta.ppf(latin_hypercube[:, idx], a, b) + latin_hypercube[:, idx] = beta.ppf( + latin_hypercube[:, idx], a, b) case "Gamma": shape, scale = params latin_hypercube[:, idx] = gamma.ppf( @@ -301,7 +305,8 @@ def validate_parameters( """ valid_strategy = ["chaospy", "scipy", "pydoe2"] - valid_distribution = ["Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma"] + valid_distribution = ["Uniform", "Normal", + "LogNormal", "Triangle", "Beta", "Gamma"] # verifying samples and distribution_params if samples is None: @@ -320,7 +325,8 @@ def validate_parameters( param = value.get("args") if dist_type is None or len(param) == 0: - raise ValueError(f"assign a list of parameters to distribution_params") + raise ValueError( + f"assign a list of parameters to distribution_params") if dist_type not in valid_distribution: raise ValueError( @@ -340,8 +346,9 @@ def validate_parameters( ) if dist_type in ["Normal", "Uniform", "Beta", "Gamma"] and len(param) != 2: - raise ValueError(f"{dist_type} distribution must have 2 parameters") - elif dist_type == "LogNormal" and len(dist_type) == 1: + raise ValueError( + f"{dist_type} distribution must have 2 parameters") + elif dist_type == "LogNormal" and len(param) != 1: raise ValueError(f"{dist_type} must have a single parameter") # handling having 0 as values in Beta and Gamma From d5dcb4e2cef317a9a0a6ecb449fd87bb99eacc42 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 22 Dec 2023 09:15:09 +0100 Subject: [PATCH 15/35] fixed lognormal shape --- scripts/monte_carlo.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 6ddd21b09..c953fb9da 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -255,20 +255,16 @@ def rescale_distribution( ) case "Normal": mean, std = params - latin_hypercube[:, idx] = norm.ppf( - latin_hypercube[:, idx], mean, std) + latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) case "LogNormal": shape = params[0] - latin_hypercube[:, idx] = lognorm.ppf( - latin_hypercube[:, idx], s=shape) + latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=shape) case "Triangle": tri_mean = np.mean(params) - latin_hypercube[:, idx] = triang.ppf( - latin_hypercube[:, idx], tri_mean) + latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) case "Beta": a, b = params - latin_hypercube[:, idx] = beta.ppf( - latin_hypercube[:, idx], a, b) + latin_hypercube[:, idx] = beta.ppf(latin_hypercube[:, idx], a, b) case "Gamma": shape, scale = params latin_hypercube[:, idx] = gamma.ppf( @@ -305,8 +301,7 @@ def validate_parameters( """ valid_strategy = ["chaospy", "scipy", "pydoe2"] - valid_distribution = ["Uniform", "Normal", - "LogNormal", "Triangle", "Beta", "Gamma"] + valid_distribution = ["Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma"] # verifying samples and distribution_params if samples is None: @@ -325,8 +320,7 @@ def validate_parameters( param = value.get("args") if dist_type is None or len(param) == 0: - raise ValueError( - f"assign a list of parameters to distribution_params") + raise ValueError(f"assign a list of parameters to distribution_params") if dist_type not in valid_distribution: raise ValueError( @@ -346,8 +340,7 @@ def validate_parameters( ) if dist_type in ["Normal", "Uniform", "Beta", "Gamma"] and len(param) != 2: - raise ValueError( - f"{dist_type} distribution must have 2 parameters") + raise ValueError(f"{dist_type} distribution must have 2 parameters") elif dist_type == "LogNormal" and len(param) != 1: raise ValueError(f"{dist_type} must have a single parameter") From 9090fc3121000818e07018bd52f9620f36f151eb Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 22 Dec 2023 09:40:32 +0100 Subject: [PATCH 16/35] reset monte-carlo option to falss --- config.default.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config.default.yaml b/config.default.yaml index 46b0d1d1e..2589ac25d 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -359,7 +359,7 @@ monte_carlo: # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object options: - add_to_snakefile: true + add_to_snakefile: false samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty From b73c295d9ae1011455a55624f816f2bae490a649 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Wed, 3 Jan 2024 19:21:41 +0100 Subject: [PATCH 17/35] fixed case to lowerletter --- config.default.yaml | 34 +++++++++--------- config.tutorial.yaml | 34 +++++++++--------- scripts/monte_carlo.py | 81 +++++++++++++++++++++++------------------- 3 files changed, 78 insertions(+), 71 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 444f39404..c329b5a73 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -340,18 +340,18 @@ costs: monte_carlo: - # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html - # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html - # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html - # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html - # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html - # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - # [mean, std] for Normal - # [shape] for LogNormal - # [lower_bound, upper_bound] for Uniform - # [lower_bound, midpoint, upper_bound] for Triangle - # [alpha, beta] for Beta - # [shape, scale] for Gamma + # uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html + # normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html + # lognormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html + # triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html + # beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html + # gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html + # [mean, std] for normal + # [shape] for lognormal + # [lower_bound, upper_bound] for uniform + # [lower_bound, midpoint, upper_bound] for triangle + # [alpha, beta] for beta + # [shape, scale] for gamma # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max(): [-3000MW 0MW]"" # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] # ... user can add flexibly more features for the Monte-Carlo sampling @@ -364,13 +364,13 @@ monte_carlo: seed: 42 # set seedling for reproducibilty uncertainties: loads_t.p_set: - type: Uniform + type: uniform args: [0, 1] - generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: - type: LogNormal + generators_t.p_max_pu.loc[:, n.generators.carrier == "onwind"]: + type: lognormal args: [1.5] - generators_t.p_min_pu.loc[:, n.generators.carrier == "wind"]: - type: Beta + generators_t.p_min_pu.loc[:, n.generators.carrier == "onwind"]: + type: beta args: [0.5, 2] diff --git a/config.tutorial.yaml b/config.tutorial.yaml index a9323e173..cd482610f 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -336,18 +336,18 @@ costs: monte_carlo: - # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html - # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html - # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html - # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html - # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html - # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - # [mean, std] for Normal - # [shape] for LogNormal - # [lower_bound, upper_bound] for Uniform - # [lower_bound, midpoint, upper_bound] for Triangle - # [alpha, beta] for Beta - # [shape, scale] for Gamma + # uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html + # normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html + # lognormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html + # triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html + # beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html + # gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html + # [mean, std] for normal + # [shape] for lognormal + # [lower_bound, upper_bound] for uniform + # [lower_bound, midpoint, upper_bound] for triangle + # [alpha, beta] for beta + # [shape, scale] for gamma # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max(): [-3000MW 0MW]"" # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] # ... user can add flexibly more features for the Monte-Carlo sampling @@ -360,13 +360,13 @@ monte_carlo: seed: 42 # set seedling for reproducibilty uncertainties: loads_t.p_set: - type: Uniform + type: uniform args: [0, 1] - generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: - type: LogNormal + generators_t.p_max_pu.loc[:, n.generators.carrier == "onwind"]: + type: lognormal args: [1.5] - generators_t.p_min_pu.loc[:, n.generators.carrier == "wind"]: - type: Beta + generators_t.p_min_pu.loc[:, n.generators.carrier == "onwind"]: + type: beta args: [0.5, 2] diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index c953fb9da..fc8c389f9 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -15,18 +15,18 @@ monte_carlo: options: add_to_snakefile: false - # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html - # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html - # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html - # Triangle: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html - # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html - # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - distribution: "Uniform" # "Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma" - # [mean, std] for Normal and LogNormal - # [lower_bound, upper_bound] for Uniform - # [lower_bound, midpoint, upper_bound] for Triangle - # [alpha, beta] for Beta - # [shape, scale] for Gamma + # uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html + # normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html + # lognormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html + # triangle: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html + # beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html + # gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html + distribution: "uniform" # "uniform", "normal", "lognormal", "triangle", "beta", "gamma" + # [mean, std] for normal and lognormal + # [lower_bound, upper_bound] for uniform + # [lower_bound, midpoint, upper_bound] for triangle + # [alpha, beta] for beta + # [shape, scale] for gamma distribution_params: [0,1] samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported @@ -229,15 +229,15 @@ def rescale_distribution( - np.array: Rescaled Latin hypercube sampling with values in the range [0, 1]. Supported Distributions: - - "Uniform": No rescaling applied. - - "Normal": Rescaled using the inverse of the normal distribution function with specified mean and std. - - "LogNormal": Rescaled using the inverse of the log-normal distribution function with specified mean and std. - - "Triangle": Rescaled using the inverse of the triangular distribution function with mean calculated from given parameters. - - "Beta": Rescaled using the inverse of the beta distribution function with specified shape parameters. - - "Gamma": Rescaled using the inverse of the gamma distribution function with specified shape and scale parameters. + - "uniform": No rescaling applied. + - "normal": Rescaled using the inverse of the normal distribution function with specified mean and std. + - "lognormal": Rescaled using the inverse of the log-normal distribution function with specified mean and std. + - "triangle": Rescaled using the inverse of the triangular distribution function with mean calculated from given parameters. + - "beta": Rescaled using the inverse of the beta distribution function with specified shape parameters. + - "gamma": Rescaled using the inverse of the gamma distribution function with specified shape and scale parameters. Note: - - The function supports rescaling for Uniform, Normal, LogNormal, Triangle, Beta, and Gamma distributions. + - The function supports rescaling for uniform, normal, lognormal, triangle, beta, and gamma distributions. - The rescaled samples will have values in the range [0, 1]. """ from scipy.stats import beta, gamma, lognorm, norm, qmc, triang @@ -248,24 +248,28 @@ def rescale_distribution( params = value.get("args") match dist: - case "Uniform": + case "uniform": l_bounds, u_bounds = params latin_hypercube[:, idx] = minmax_scale( latin_hypercube[:, idx], feature_range=(l_bounds, u_bounds) ) - case "Normal": + case "normal": mean, std = params - latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) - case "LogNormal": + latin_hypercube[:, idx] = norm.ppf( + latin_hypercube[:, idx], mean, std) + case "lognormal": shape = params[0] - latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=shape) - case "Triangle": + latin_hypercube[:, idx] = lognorm.ppf( + latin_hypercube[:, idx], s=shape) + case "triangle": tri_mean = np.mean(params) - latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) - case "Beta": + latin_hypercube[:, idx] = triang.ppf( + latin_hypercube[:, idx], tri_mean) + case "beta": a, b = params - latin_hypercube[:, idx] = beta.ppf(latin_hypercube[:, idx], a, b) - case "Gamma": + latin_hypercube[:, idx] = beta.ppf( + latin_hypercube[:, idx], a, b) + case "gamma": shape, scale = params latin_hypercube[:, idx] = gamma.ppf( latin_hypercube[:, idx], shape, scale @@ -301,7 +305,8 @@ def validate_parameters( """ valid_strategy = ["chaospy", "scipy", "pydoe2"] - valid_distribution = ["Uniform", "Normal", "LogNormal", "Triangle", "Beta", "Gamma"] + valid_distribution = ["uniform", "normal", + "lognormal", "triangle", "beta", "gamma"] # verifying samples and distribution_params if samples is None: @@ -320,14 +325,15 @@ def validate_parameters( param = value.get("args") if dist_type is None or len(param) == 0: - raise ValueError(f"assign a list of parameters to distribution_params") + raise ValueError( + f"assign a list of parameters to distribution_params") if dist_type not in valid_distribution: raise ValueError( f"Unsupported Distribution : {dist_type}. Choose from {valid_distribution}" ) - if dist_type == "Triangle": + if dist_type == "triangle": if len(param) == 2: print( f"{dist_type} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" @@ -339,13 +345,14 @@ def validate_parameters( f"{dist_type} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" ) - if dist_type in ["Normal", "Uniform", "Beta", "Gamma"] and len(param) != 2: - raise ValueError(f"{dist_type} distribution must have 2 parameters") - elif dist_type == "LogNormal" and len(param) != 1: + if dist_type in ["normal", "uniform", "beta", "gamma"] and len(param) != 2: + raise ValueError( + f"{dist_type} distribution must have 2 parameters") + elif dist_type == "lognormal" and len(param) != 1: raise ValueError(f"{dist_type} must have a single parameter") - # handling having 0 as values in Beta and Gamma - if dist_type in ["Beta", "Gamma"]: + # handling having 0 as values in beta and gamma + if dist_type in ["beta", "gamma"]: if np.min(param) <= 0: raise ValueError( f"{dist_type} distribution cannot have values lower than zero in parameters" From 7af5c4e803b4d357e6d2a64c0fbade55a249300c Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Wed, 3 Jan 2024 19:22:15 +0100 Subject: [PATCH 18/35] fixed case to lowerletter --- scripts/monte_carlo.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index fc8c389f9..2ad6e5ded 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -255,20 +255,16 @@ def rescale_distribution( ) case "normal": mean, std = params - latin_hypercube[:, idx] = norm.ppf( - latin_hypercube[:, idx], mean, std) + latin_hypercube[:, idx] = norm.ppf(latin_hypercube[:, idx], mean, std) case "lognormal": shape = params[0] - latin_hypercube[:, idx] = lognorm.ppf( - latin_hypercube[:, idx], s=shape) + latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=shape) case "triangle": tri_mean = np.mean(params) - latin_hypercube[:, idx] = triang.ppf( - latin_hypercube[:, idx], tri_mean) + latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) case "beta": a, b = params - latin_hypercube[:, idx] = beta.ppf( - latin_hypercube[:, idx], a, b) + latin_hypercube[:, idx] = beta.ppf(latin_hypercube[:, idx], a, b) case "gamma": shape, scale = params latin_hypercube[:, idx] = gamma.ppf( @@ -305,8 +301,7 @@ def validate_parameters( """ valid_strategy = ["chaospy", "scipy", "pydoe2"] - valid_distribution = ["uniform", "normal", - "lognormal", "triangle", "beta", "gamma"] + valid_distribution = ["uniform", "normal", "lognormal", "triangle", "beta", "gamma"] # verifying samples and distribution_params if samples is None: @@ -325,8 +320,7 @@ def validate_parameters( param = value.get("args") if dist_type is None or len(param) == 0: - raise ValueError( - f"assign a list of parameters to distribution_params") + raise ValueError(f"assign a list of parameters to distribution_params") if dist_type not in valid_distribution: raise ValueError( @@ -346,8 +340,7 @@ def validate_parameters( ) if dist_type in ["normal", "uniform", "beta", "gamma"] and len(param) != 2: - raise ValueError( - f"{dist_type} distribution must have 2 parameters") + raise ValueError(f"{dist_type} distribution must have 2 parameters") elif dist_type == "lognormal" and len(param) != 1: raise ValueError(f"{dist_type} must have a single parameter") From c2235924ecffc27990aa6bdddb59cc1ffac7e65f Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 9 Jan 2024 12:06:12 +0100 Subject: [PATCH 19/35] fixed plots --- scripts/monte_carlo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 2ad6e5ded..584adb51c 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -93,7 +93,7 @@ from solve_network import * logger = create_logger(__name__) -sns.set() +sns.set(style="whitegrid") def monte_carlo_sampling_pydoe2( @@ -419,7 +419,7 @@ def validate_parameters( # create plot for the rescaled distributions for idx in range(N_FEATURES): - sns.displot(lh[:, idx]).set( + sns.histplot(lh[:, idx], kde=True).set( title=f"{MONTE_CARLO_PYPSA_FEATURES[idx]}" ).figure.savefig(f"{MONTE_CARLO_PYPSA_FEATURES[idx]}.png", bbox_inches="tight") From c63720615691f34ca5c6ff4d988bee0c8fb58145 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 9 Jan 2024 12:44:24 +0100 Subject: [PATCH 20/35] fixed plots back to displot --- scripts/monte_carlo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 584adb51c..b22711a0b 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -419,7 +419,7 @@ def validate_parameters( # create plot for the rescaled distributions for idx in range(N_FEATURES): - sns.histplot(lh[:, idx], kde=True).set( + sns.displot(lh[:, idx], kde=True).set( title=f"{MONTE_CARLO_PYPSA_FEATURES[idx]}" ).figure.savefig(f"{MONTE_CARLO_PYPSA_FEATURES[idx]}.png", bbox_inches="tight") From 9772850fe3ae38a8a26d6471447841528f1b71b8 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 12 Jan 2024 11:42:00 +0100 Subject: [PATCH 21/35] modified clean rule --- Snakefile | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/Snakefile b/Snakefile index ebc220988..9948d364d 100644 --- a/Snakefile +++ b/Snakefile @@ -65,20 +65,12 @@ if config["custom_rules"] is not []: include: rule -rule clean_monte: - run: - try: - shell("snakemake -j 1 solve_all_networks_monte --delete-all-output") - except: - pass - shell("snakemake -j 1 run_all_scenarios --delete-all-output") - - rule clean: run: try: shell("snakemake -j 1 solve_all_networks --delete-all-output") except: + shell("snakemake -j 1 solve_all_networks_monte --delete-all-output") pass shell("snakemake -j 1 run_all_scenarios --delete-all-output") From 5f0f74dc3fc94d026825d1699afa09a0999af091 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 12 Jan 2024 11:44:51 +0100 Subject: [PATCH 22/35] revert defaults --- config.default.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config.default.yaml b/config.default.yaml index c329b5a73..335a86592 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -25,7 +25,7 @@ enable: custom_rules: [] # Default empty [] or link to custom rule file e.g. ["my_folder/my_rules.smk"] that add rules to Snakefile run: - name: "NG_BJ" # use this to keep track of runs with different settings + name: "" # use this to keep track of runs with different settings shared_cutouts: true # set to true to share the default cutout(s) across runs # Note: value false requires build_cutout to be enabled From 34346bef367430d0551af0881dc80dee9e1f8589 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 12 Jan 2024 11:47:51 +0100 Subject: [PATCH 23/35] added comments to enable monte-carlo --- config.default.yaml | 2 +- config.tutorial.yaml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 335a86592..ea0187e00 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -358,7 +358,7 @@ monte_carlo: # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object options: - add_to_snakefile: false + add_to_snakefile: false # when set to true, monte-carlo is enabled samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty diff --git a/config.tutorial.yaml b/config.tutorial.yaml index cd482610f..705039b2e 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -354,7 +354,7 @@ monte_carlo: # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object options: - add_to_snakefile: false + add_to_snakefile: false # when set to true, monte-carlo is enabled samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty From bb31bfe112e1e9581186dab55bda187eed286796 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 12 Jan 2024 11:53:44 +0100 Subject: [PATCH 24/35] commented out plots for rescaled distribution --- scripts/monte_carlo.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index b22711a0b..d3ed1fe2f 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -417,11 +417,11 @@ def validate_parameters( N_FEATURES, SAMPLES, UNCERTAINTIES_VALUES, seed=SEED, rule="latin_hypercube" ) - # create plot for the rescaled distributions - for idx in range(N_FEATURES): - sns.displot(lh[:, idx], kde=True).set( - title=f"{MONTE_CARLO_PYPSA_FEATURES[idx]}" - ).figure.savefig(f"{MONTE_CARLO_PYPSA_FEATURES[idx]}.png", bbox_inches="tight") + # create plot for the rescaled distributions (for development usage, commented by default) + # for idx in range(N_FEATURES): + # sns.displot(lh[:, idx], kde=True).set( + # title=f"{MONTE_CARLO_PYPSA_FEATURES[idx]}" + # ).figure.savefig(f"{MONTE_CARLO_PYPSA_FEATURES[idx]}.png", bbox_inches="tight") # MONTE-CARLO MODIFICATIONS ### From 127d797565cd46242df9f6eafcfb67bf88d2df4a Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Mon, 15 Jan 2024 08:50:45 +0100 Subject: [PATCH 25/35] updated comments for MC --- config.default.yaml | 32 ++++++++++++++------------------ config.tutorial.yaml | 30 +++++++++++------------------- 2 files changed, 25 insertions(+), 37 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index ea0187e00..0038af4ad 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -340,28 +340,20 @@ costs: monte_carlo: - # uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html - # normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html - # lognormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html - # triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html - # beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html - # gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - # [mean, std] for normal - # [shape] for lognormal - # [lower_bound, upper_bound] for uniform - # [lower_bound, midpoint, upper_bound] for triangle - # [alpha, beta] for beta - # [shape, scale] for gamma - # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max(): [-3000MW 0MW]"" - # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] - # ... user can add flexibly more features for the Monte-Carlo sampling - # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! - # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object + # Description: Specify Monte Carlo sampling options for uncertainty analysis. + # Define the option list for Monte Carlo sampling. + # Make sure add_to_snakefile is set to true to enable Monte-Carlo options: - add_to_snakefile: false # when set to true, monte-carlo is enabled + add_to_snakefile: false # When set to true, enables Monte Carlo sampling samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty + + # Different distributions can be specified for various PyPSA network parameters. + # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. + # More info on the distributions are documented in the Chaospy reference guide... + # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html + # ... users can add flexibly more features for the Monte-Carlo sampling uncertainties: loads_t.p_set: type: uniform @@ -373,6 +365,10 @@ monte_carlo: type: beta args: [0.5, 2] + # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max()": [-3000MW 0MW] + # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] + + solving: options: diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 705039b2e..1fa46dcec 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -336,28 +336,19 @@ costs: monte_carlo: - # uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html - # normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html - # lognormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html - # triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html - # beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html - # gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - # [mean, std] for normal - # [shape] for lognormal - # [lower_bound, upper_bound] for uniform - # [lower_bound, midpoint, upper_bound] for triangle - # [alpha, beta] for beta - # [shape, scale] for gamma - # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max(): [-3000MW 0MW]"" - # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] - # ... user can add flexibly more features for the Monte-Carlo sampling - # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! - # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object + # Description: Specify Monte Carlo sampling options for uncertainty analysis. + # Define the option list for Monte Carlo sampling. + # Make sure add_to_snakefile is set to true to enable Monte-Carlo options: - add_to_snakefile: false # when set to true, monte-carlo is enabled + add_to_snakefile: false # When set to true, enables Monte Carlo sampling samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty + # Different distributions can be specified for various PyPSA network parameters. + # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. + # More info on the distributions are documented in the Chaospy reference guide... + # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html + # ... users can add flexibly more features for the Monte-Carlo sampling uncertainties: loads_t.p_set: type: uniform @@ -368,7 +359,8 @@ monte_carlo: generators_t.p_min_pu.loc[:, n.generators.carrier == "onwind"]: type: beta args: [0.5, 2] - + # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max()": [-3000MW 0MW] + # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] solving: options: From 8079f2207e27ef1454f738ce355682d85b0f2062 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Fri, 19 Jan 2024 15:04:30 +0100 Subject: [PATCH 26/35] more comments to config documentation --- config.default.yaml | 9 +++++---- config.tutorial.yaml | 7 +++++-- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 0038af4ad..dc0da722f 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -348,12 +348,14 @@ monte_carlo: samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty - - # Different distributions can be specified for various PyPSA network parameters. + # Different distributions can be specified for various PyPSA network object. # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html - # ... users can add flexibly more features for the Monte-Carlo sampling + # ... users can add flexibly more features for the Monte-Carlo sampling using the description below + # {pypsa network object, e.g. "loads_t.p_set"}: + # type: {any distribution among: "uniform", "normal", "lognormal", "triangle", "beta" and "gamma"} + # args: {arguments passed as a list depending on the distribution, see arguments description in Chaospy reference guide} uncertainties: loads_t.p_set: type: uniform @@ -364,7 +366,6 @@ monte_carlo: generators_t.p_min_pu.loc[:, n.generators.carrier == "onwind"]: type: beta args: [0.5, 2] - # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max()": [-3000MW 0MW] # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW] diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 1fa46dcec..85b2f66ff 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -344,11 +344,14 @@ monte_carlo: samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty - # Different distributions can be specified for various PyPSA network parameters. + # Different distributions can be specified for various PyPSA network object. # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html - # ... users can add flexibly more features for the Monte-Carlo sampling + # ... users can add flexibly more features for the Monte-Carlo sampling using the description below + # {pypsa network object, e.g. "loads_t.p_set"}: + # type: {any distribution among: "uniform", "normal", "lognormal", "triangle", "beta" and "gamma"} + # args: {arguments passed as a list depending on the distribution, see arguments description in Chaospy reference guide} uncertainties: loads_t.p_set: type: uniform From 4862e9beec88854ede802d0898fd2297e401f6f3 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 23 Jan 2024 14:45:07 +0100 Subject: [PATCH 27/35] include monte-carlo into api-reference --- doc/api_reference.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/doc/api_reference.rst b/doc/api_reference.rst index b5832936b..e24cf46ce 100644 --- a/doc/api_reference.rst +++ b/doc/api_reference.rst @@ -103,3 +103,9 @@ make_statistics .. automodule:: make_statistics :members: + +monte-carlo +------------------------------- + +.. automodule:: monte-carlo + :members: From c8e5e2f3a12e941b684c7301a6bd08edc40dd3e2 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Tue, 23 Jan 2024 14:49:33 +0100 Subject: [PATCH 28/35] modified description --- scripts/monte_carlo.py | 49 +++++++++++++++++++++--------------------- 1 file changed, 25 insertions(+), 24 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index d3ed1fe2f..c3ec63330 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -13,31 +13,32 @@ .. code:: yaml monte_carlo: + # Description: Specify Monte Carlo sampling options for uncertainty analysis. + # Define the option list for Monte Carlo sampling. + # Make sure add_to_snakefile is set to true to enable Monte-Carlo options: - add_to_snakefile: false - # uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html - # normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html - # lognormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html - # triangle: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html - # beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html - # gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html - distribution: "uniform" # "uniform", "normal", "lognormal", "triangle", "beta", "gamma" - # [mean, std] for normal and lognormal - # [lower_bound, upper_bound] for uniform - # [lower_bound, midpoint, upper_bound] for triangle - # [alpha, beta] for beta - # [shape, scale] for gamma - distribution_params: [0,1] - samples: 4 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number - sampling_strategy: "scipy" # "pydoe2", "chaospy", "scipy", packages that are supported - pypsa_standard: - # User can add here flexibly more features for the Monte-Carlo sampling. - # Given as "key: value" format - # Key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering! - # Value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object - loads_t.p_set: [0.9, 1.1] - # generators_t.p_max_pu.loc[:, n.generators.carrier == "wind"]: [0.9, 1.1] - # generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]: [0.9, 1.1] + add_to_snakefile: false # When set to true, enables Monte Carlo sampling + samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number + sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported + seed: 42 # set seedling for reproducibilty + # Different distributions can be specified for various PyPSA network object. + # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. + # More info on the distributions are documented in the Chaospy reference guide... + # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html + # ... users can add flexibly more features for the Monte-Carlo sampling using the description below + # {pypsa network object, e.g. "loads_t.p_set"}: + # type: {any distribution among: "uniform", "normal", "lognormal", "triangle", "beta" and "gamma"} + # args: {arguments passed as a list depending on the distribution, see arguments description in Chaospy reference guide} + uncertainties: + loads_t.p_set: + type: uniform + args: [0, 1] + generators_t.p_max_pu.loc[:, n.generators.carrier == "onwind"]: + type: lognormal + args: [1.5] + generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]: + type: beta + args: [0.5, 2] .. seealso:: Documentation of the configuration file ``config.yaml`` at :ref:`_monte_cf` From 171bde7a5f64e9e5b3712f3feef0c6a32fb62a55 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Wed, 24 Jan 2024 11:06:13 +0100 Subject: [PATCH 29/35] fix monte-carlo --- doc/api_reference.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/api_reference.rst b/doc/api_reference.rst index e24cf46ce..5cf8fdf3a 100644 --- a/doc/api_reference.rst +++ b/doc/api_reference.rst @@ -104,8 +104,8 @@ make_statistics .. automodule:: make_statistics :members: -monte-carlo +monte_carlo ------------------------------- -.. automodule:: monte-carlo +.. automodule:: monte_carlo :members: From cafc172ae339b6d73d250d6278a3c224612836a8 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Wed, 24 Jan 2024 11:08:27 +0100 Subject: [PATCH 30/35] remove old comments from description --- scripts/monte_carlo.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index c3ec63330..799a25095 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -25,7 +25,7 @@ # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html - # ... users can add flexibly more features for the Monte-Carlo sampling using the description below + # ... users can flexibly add more features for the Monte-Carlo sampling in the uncertainties tab using the description below # {pypsa network object, e.g. "loads_t.p_set"}: # type: {any distribution among: "uniform", "normal", "lognormal", "triangle", "beta" and "gamma"} # args: {arguments passed as a list depending on the distribution, see arguments description in Chaospy reference guide} @@ -41,7 +41,7 @@ args: [0.5, 2] .. seealso:: - Documentation of the configuration file ``config.yaml`` at :ref:`_monte_cf` + Documentation of the configuration file ``config.yaml`` at :ref:`monte_cf` Inputs ------ @@ -68,13 +68,11 @@ supported by the packages pyDOE2, chaospy and scipy. This should give users the freedom to explore alternative approaches. The orthogonal latin hypercube sampling is thereby found as most performant, hence, implemented here. Sampling the multi-dimensional uncertainty space is relatively -easy. It only requires two things: The number of *samples* (e.g. PyPSA networks) and *features* (e.g. -load or solar timeseries). This results in an experimental design of the dimension (samples X features). +easy. It only requires two things: The number of *samples* (defines the number of total networks to +be optimized) and *features* (pypsa network object e.g loads_t.p_set or generators_t.p_max_pu). This +results in an experimental design of the dimension (samples X features). -Additionally, upper and lower bounds *per feature* need to be provided such that the experimental -design can be scaled accordingly. Currently the user can define uncertainty ranges e.g. bounds, -for all PyPSA objects that are `int` or `float`. Boolean values could be used but require testing. -The experimental design `lhs_scaled` (dimension: sample X features) is then used to modify the PyPSA +The experimental design `lh` (dimension: sample X features) is used to modify the PyPSA networks. Thereby, this script creates samples x amount of networks. The iterators comes from the wildcard {unc}, which is described in the config.yaml and created in the Snakefile as a range from 0 to (total number of) SAMPLES. From 7caed111223166a02452b960338704c1bf2070b7 Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Wed, 24 Jan 2024 12:17:35 +0100 Subject: [PATCH 31/35] adjusted comments for readthedocs --- config.default.yaml | 5 ++++ config.tutorial.yaml | 5 ++++ scripts/monte_carlo.py | 52 +++++++++++++++++++++++------------------- 3 files changed, 39 insertions(+), 23 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index dc0da722f..e6bb890cd 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -350,6 +350,11 @@ monte_carlo: seed: 42 # set seedling for reproducibilty # Different distributions can be specified for various PyPSA network object. # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. + # [mean, std] for Normal and LogNormal + # [lower_bound, upper_bound] for Uniform + # [mid_point (between 0 - 1)] for Triangle + # [alpha, beta] for Beta + # [shape, scale] for Gamma # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html # ... users can add flexibly more features for the Monte-Carlo sampling using the description below diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 85b2f66ff..b99d84d25 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -346,6 +346,11 @@ monte_carlo: seed: 42 # set seedling for reproducibilty # Different distributions can be specified for various PyPSA network object. # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. + # [mean, std] for Normal and LogNormal + # [lower_bound, upper_bound] for Uniform + # [mid_point (between 0 - 1)] for Triangle + # [alpha, beta] for Beta + # [shape, scale] for Gamma # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html # ... users can add flexibly more features for the Monte-Carlo sampling using the description below diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 799a25095..2fbf263d8 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -23,6 +23,11 @@ seed: 42 # set seedling for reproducibilty # Different distributions can be specified for various PyPSA network object. # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. + # [mean, std] for Normal and LogNormal + # [lower_bound, upper_bound] for Uniform + # [mid_point (between 0 - 1)] for Triangle + # [alpha, beta] for Beta + # [shape, scale] for Gamma # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html # ... users can flexibly add more features for the Monte-Carlo sampling in the uncertainties tab using the description below @@ -179,9 +184,9 @@ def monte_carlo_sampling_scipy( options: - Center the point within the multi-dimensional grid, centered=True - - optimization scheme, optimization="random-cd" - - strength=1, classical LHS - - strength=2, performant orthogonal LHS, requires the sample to be a prime or square of a prime e.g. sqr(121)=11 + - Optimization scheme, optimization="random-cd" + - Strength=1, classical LHS + - Strength=2, performant orthogonal LHS, requires the sample to be square of a prime e.g. sq(11)=121 Options could be combined to produce an optimized centered orthogonal array based LHS. After optimization, the result would not be guaranteed to be of strength 2. @@ -216,26 +221,31 @@ def rescale_distribution( ) -> np.ndarray: """ Rescales a Latin hypercube sampling (LHS) using specified distribution - parameters. + parameters. More information on the distributions can be found + here https://docs.scipy.org/doc/scipy/reference/stats.html + + **Parameters**: - Parameters: - latin_hypercube (np.array): The Latin hypercube sampling to be rescaled. - uncertainties_values (list): List of dictionaries containing distribution information. - Each dictionary should have 'type' key specifying the distribution type and 'args' key - containing parameters specific to the chosen distribution. + Each dictionary should have 'type' key specifying the distribution type and 'args' key + containing parameters specific to the chosen distribution. + + **Returns**: - Returns: - np.array: Rescaled Latin hypercube sampling with values in the range [0, 1]. - Supported Distributions: - - "uniform": No rescaling applied. + **Supported Distributions**: + + - "uniform": Rescaled to the specified lower and upper bounds. - "normal": Rescaled using the inverse of the normal distribution function with specified mean and std. - "lognormal": Rescaled using the inverse of the log-normal distribution function with specified mean and std. - "triangle": Rescaled using the inverse of the triangular distribution function with mean calculated from given parameters. - "beta": Rescaled using the inverse of the beta distribution function with specified shape parameters. - "gamma": Rescaled using the inverse of the gamma distribution function with specified shape and scale parameters. - Note: + **Note**: + - The function supports rescaling for uniform, normal, lognormal, triangle, beta, and gamma distributions. - The rescaled samples will have values in the range [0, 1]. """ @@ -259,8 +269,8 @@ def rescale_distribution( shape = params[0] latin_hypercube[:, idx] = lognorm.ppf(latin_hypercube[:, idx], s=shape) case "triangle": - tri_mean = np.mean(params) - latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], tri_mean) + mid_point = params[0] + latin_hypercube[:, idx] = triang.ppf(latin_hypercube[:, idx], mid_point) case "beta": a, b = params latin_hypercube[:, idx] = beta.ppf(latin_hypercube[:, idx], a, b) @@ -285,7 +295,8 @@ def validate_parameters( user through the config file needs to be validated before proceeding to perform monte-carlo simulations. - Parameters: + **Parameters**: + - sampling_strategy: str The chosen sampling strategy from chaospy, scipy and pydoe2 - samples: int @@ -295,7 +306,8 @@ def validate_parameters( - distribution_params: list The parameters associated with the probability distribution. - Raises: + **Raises**: + - ValueError: If the parameters are invalid for the specified distribution. """ @@ -327,15 +339,9 @@ def validate_parameters( ) if dist_type == "triangle": - if len(param) == 2: - print( - f"{dist_type} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" - ) - # use the mean as the middle value - param.insert(1, np.mean(param)) - elif len(param) != 3: + if len(param) != 1: raise ValueError( - f"{dist_type} distribution has to be 3 parameters in the order of [lower_bound, mid_range, upper_bound]" + f"{dist_type} distribution has to be 1 parameter [midpoint_between_0_and_1]" ) if dist_type in ["normal", "uniform", "beta", "gamma"] and len(param) != 2: From d0fa46935635f7ad374ca21f81c324bd8d77ca6b Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Wed, 24 Jan 2024 13:29:50 +0100 Subject: [PATCH 32/35] fixed tables for monte-carlo in docs --- config.default.yaml | 9 +++------ config.tutorial.yaml | 9 +++------ doc/configtables/monte-carlo.csv | 15 +++++++++------ scripts/monte_carlo.py | 7 ++----- 4 files changed, 17 insertions(+), 23 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index e6bb890cd..b094b1896 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -350,11 +350,8 @@ monte_carlo: seed: 42 # set seedling for reproducibilty # Different distributions can be specified for various PyPSA network object. # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. - # [mean, std] for Normal and LogNormal - # [lower_bound, upper_bound] for Uniform - # [mid_point (between 0 - 1)] for Triangle - # [alpha, beta] for Beta - # [shape, scale] for Gamma + # normal: [mean, std], lognormal: [mean, std], uniform: [lower_bound, upper_bound], + # triangle: [mid_point (between 0 - 1)], beta: [alpha, beta], gamma: [shape, scale] # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html # ... users can add flexibly more features for the Monte-Carlo sampling using the description below @@ -368,7 +365,7 @@ monte_carlo: generators_t.p_max_pu.loc[:, n.generators.carrier == "onwind"]: type: lognormal args: [1.5] - generators_t.p_min_pu.loc[:, n.generators.carrier == "onwind"]: + generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]: type: beta args: [0.5, 2] # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max()": [-3000MW 0MW] diff --git a/config.tutorial.yaml b/config.tutorial.yaml index b99d84d25..d3ca6b96d 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -346,11 +346,8 @@ monte_carlo: seed: 42 # set seedling for reproducibilty # Different distributions can be specified for various PyPSA network object. # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. - # [mean, std] for Normal and LogNormal - # [lower_bound, upper_bound] for Uniform - # [mid_point (between 0 - 1)] for Triangle - # [alpha, beta] for Beta - # [shape, scale] for Gamma + # normal: [mean, std], lognormal: [mean, std], uniform: [lower_bound, upper_bound], + # triangle: [mid_point (between 0 - 1)], beta: [alpha, beta], gamma: [shape, scale] # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html # ... users can add flexibly more features for the Monte-Carlo sampling using the description below @@ -364,7 +361,7 @@ monte_carlo: generators_t.p_max_pu.loc[:, n.generators.carrier == "onwind"]: type: lognormal args: [1.5] - generators_t.p_min_pu.loc[:, n.generators.carrier == "onwind"]: + generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]: type: beta args: [0.5, 2] # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max()": [-3000MW 0MW] diff --git a/doc/configtables/monte-carlo.csv b/doc/configtables/monte-carlo.csv index efde31bd0..bc922f2f3 100644 --- a/doc/configtables/monte-carlo.csv +++ b/doc/configtables/monte-carlo.csv @@ -1,7 +1,10 @@ ,Unit,Values,Description -options,,, --- add_to_snakemake,,true or false,Add rule to Snakefile or not --- samples,,"int** read description!", "Defines the number of total sample networks that will be later optimized. if orthogonal latin hypercube sampling is applied (default option in scripts), then a prime number, that is not a quadrat of a prime, needs to be chosen. E.g. 7 not 9 (3^2)" --- sampling_strategy,,"Any subset of {""pydoe2"", ""chaospy"", ""scipy""}",Current supported packages to create an experimental design -pypsa_standard,,, --- ,MW/MWh,"[l_bound, u_bound] or empty []","`Key` is a dynamic PyPSA object that allows to access any pypsa object such as `loads_t.p_set` or the max. wind generation per hour `generators_t.p_max_pu.loc[:, n.generators.carrier == ""wind""]`. `Values` or bounds are multiplication for each object." +**options**,,, +add_to_snakemake,,"true or false","Set to true to enable Monte-Carlo" +samples,,"int","Defines the number of total sample networks that will be optimized. If the chosen sampling strategy is scipy, then a square of a prime number needs to be chosen. E.g. 49 which is (7^2)" +sampling_strategy,,"Any subset of {pydoe2, chaospy, scipy}","Current supported packages to create an experimental design" +seed,,"int","Allows experimentation to be reproduced easily" +**uncertainties**,,, +,MW/MWh,,"`Key` is a dynamic PyPSA object that allows to access any pypsa object such as `loads_t.p_set` or the max. wind generation per hour `generators_t.p_max_pu.loc[:, n.generators.carrier == ""wind""]`. `Values` or bounds are multiplication for each object." +type,,"str","Defines the distribution for the chosen pypsa.object parameter. Distribution can be either uniform, normal, lognormal, triangle, beta or gamma" +args,,"list","Defines parameters for the chosen distribution. [mean, std] for normal and lognormal, [lower_bound, upper_bound] for uniform, [mid_point (between 0 - 1)] for triangle, [alpha, beta] for beta, [shape, scale] for gamma" diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 2fbf263d8..3f6fc6243 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -23,11 +23,8 @@ seed: 42 # set seedling for reproducibilty # Different distributions can be specified for various PyPSA network object. # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. - # [mean, std] for Normal and LogNormal - # [lower_bound, upper_bound] for Uniform - # [mid_point (between 0 - 1)] for Triangle - # [alpha, beta] for Beta - # [shape, scale] for Gamma + # normal: [mean, std], lognormal: [mean, std], uniform: [lower_bound, upper_bound], + # triangle: [mid_point (between 0 - 1)], beta: [alpha, beta], gamma: [shape, scale] # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html # ... users can flexibly add more features for the Monte-Carlo sampling in the uncertainties tab using the description below From f4053c1262ce6d1eadd67d683fe465bc3b2b2323 Mon Sep 17 00:00:00 2001 From: Emmanuel Bolarinwa Date: Wed, 24 Jan 2024 21:41:23 +0100 Subject: [PATCH 33/35] Update config.default.yaml Co-authored-by: Davide Fioriti <67809479+davide-f@users.noreply.github.com> --- config.default.yaml | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index b094b1896..d70c3e88f 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -348,16 +348,18 @@ monte_carlo: samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty - # Different distributions can be specified for various PyPSA network object. - # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. + # Uncertanties on any PyPSA object are specified by declaring the specific PyPSA object under the key 'uncertainties'. + # For each PyPSA object, the 'type' and 'args' keys represent the type of distribution and its argument, respectively. + # Supported distributions types are uniform, normal, lognormal, triangle, beta and gamma. + # The arguments of the distribution are passed using the key 'args' as follows, tailored by distribution type # normal: [mean, std], lognormal: [mean, std], uniform: [lower_bound, upper_bound], # triangle: [mid_point (between 0 - 1)], beta: [alpha, beta], gamma: [shape, scale] # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html - # ... users can add flexibly more features for the Monte-Carlo sampling using the description below + # An abstract example is as follows: # {pypsa network object, e.g. "loads_t.p_set"}: - # type: {any distribution among: "uniform", "normal", "lognormal", "triangle", "beta" and "gamma"} - # args: {arguments passed as a list depending on the distribution, see arguments description in Chaospy reference guide} + # type: {any supported distribution among the previous: "uniform", "normal", ...} + # args: {arguments passed as a list depending on the distribution, see the above and more at https://pypsa.readthedocs.io/} uncertainties: loads_t.p_set: type: uniform From 5741ee77ced2e0798eb8bc1629a88fb14781ca4d Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Thu, 25 Jan 2024 03:33:17 +0100 Subject: [PATCH 34/35] removed unnecessary comments --- config.tutorial.yaml | 12 +++++++----- scripts/monte_carlo.py | 13 ------------- 2 files changed, 7 insertions(+), 18 deletions(-) diff --git a/config.tutorial.yaml b/config.tutorial.yaml index d3ca6b96d..8ff9630cb 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -344,16 +344,18 @@ monte_carlo: samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty - # Different distributions can be specified for various PyPSA network object. - # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. + # Uncertanties on any PyPSA object are specified by declaring the specific PyPSA object under the key 'uncertainties'. + # For each PyPSA object, the 'type' and 'args' keys represent the type of distribution and its argument, respectively. + # Supported distributions types are uniform, normal, lognormal, triangle, beta and gamma. + # The arguments of the distribution are passed using the key 'args' as follows, tailored by distribution type # normal: [mean, std], lognormal: [mean, std], uniform: [lower_bound, upper_bound], # triangle: [mid_point (between 0 - 1)], beta: [alpha, beta], gamma: [shape, scale] # More info on the distributions are documented in the Chaospy reference guide... # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html - # ... users can add flexibly more features for the Monte-Carlo sampling using the description below + # An abstract example is as follows: # {pypsa network object, e.g. "loads_t.p_set"}: - # type: {any distribution among: "uniform", "normal", "lognormal", "triangle", "beta" and "gamma"} - # args: {arguments passed as a list depending on the distribution, see arguments description in Chaospy reference guide} + # type: {any supported distribution among the previous: "uniform", "normal", ...} + # args: {arguments passed as a list depending on the distribution, see the above and more at https://pypsa.readthedocs.io/} uncertainties: loads_t.p_set: type: uniform diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 3f6fc6243..02cafa6c1 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -13,24 +13,11 @@ .. code:: yaml monte_carlo: - # Description: Specify Monte Carlo sampling options for uncertainty analysis. - # Define the option list for Monte Carlo sampling. - # Make sure add_to_snakefile is set to true to enable Monte-Carlo options: add_to_snakefile: false # When set to true, enables Monte Carlo sampling samples: 9 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number sampling_strategy: "chaospy" # "pydoe2", "chaospy", "scipy", packages that are supported seed: 42 # set seedling for reproducibilty - # Different distributions can be specified for various PyPSA network object. - # Supported distributions for uncertainties are uniform, normal, lognormal, triangle, beta and gamma. - # normal: [mean, std], lognormal: [mean, std], uniform: [lower_bound, upper_bound], - # triangle: [mid_point (between 0 - 1)], beta: [alpha, beta], gamma: [shape, scale] - # More info on the distributions are documented in the Chaospy reference guide... - # https://chaospy.readthedocs.io/en/master/reference/distribution/index.html - # ... users can flexibly add more features for the Monte-Carlo sampling in the uncertainties tab using the description below - # {pypsa network object, e.g. "loads_t.p_set"}: - # type: {any distribution among: "uniform", "normal", "lognormal", "triangle", "beta" and "gamma"} - # args: {arguments passed as a list depending on the distribution, see arguments description in Chaospy reference guide} uncertainties: loads_t.p_set: type: uniform From 5a1b2f111490e74af6f0806d645ff88c0fc61b2f Mon Sep 17 00:00:00 2001 From: GbotemiB Date: Thu, 25 Jan 2024 16:13:01 +0100 Subject: [PATCH 35/35] fixed distribution issue --- scripts/monte_carlo.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/scripts/monte_carlo.py b/scripts/monte_carlo.py index 02cafa6c1..2338e892c 100644 --- a/scripts/monte_carlo.py +++ b/scripts/monte_carlo.py @@ -159,7 +159,6 @@ def monte_carlo_sampling_scipy( SAMPLES: int, uncertainties_values: dict, seed: int, - centered: bool = False, strength: int = 2, optimization: str = None, ) -> np.ndarray: @@ -183,7 +182,6 @@ def monte_carlo_sampling_scipy( sampler = qmc.LatinHypercube( d=N_FEATURES, - centered=centered, strength=strength, optimization=optimization, seed=seed, @@ -327,6 +325,10 @@ def validate_parameters( raise ValueError( f"{dist_type} distribution has to be 1 parameter [midpoint_between_0_and_1]" ) + elif param[0] < 0 or param[0] > 1: + raise ValueError( + f"{dist_type} distribution has to have a midpoint between 0 and 1" + ) if dist_type in ["normal", "uniform", "beta", "gamma"] and len(param) != 2: raise ValueError(f"{dist_type} distribution must have 2 parameters") @@ -397,7 +399,6 @@ def validate_parameters( SAMPLES, UNCERTAINTIES_VALUES, seed=SEED, - centered=False, strength=2, optimization=None, )