From dcc8da01edebc2c3f836678f69cbd33b4414e132 Mon Sep 17 00:00:00 2001 From: Sylvain Leclerc Date: Sun, 28 Jul 2024 14:38:45 +0200 Subject: [PATCH] Fix inconsistency with cpp impl, improve design (#3) - fix: RNG is skipped when volatility is 0 - RNG as an argument to higher level generators - Separate computation of available units and available power - add unit tests Signed-off-by: Sylvain Leclerc --- requirements.txt | 1 + src/cluster_import.py | 15 ++- src/cluster_resolve.py | 7 +- src/duration_generator.py | 115 +++++++++++++++++++ src/mersenne_twister.py | 9 +- src/random_generator.py | 48 ++++++++ src/ts_generator.py | 169 ++++++++------------------- tests/conftest.py | 7 ++ tests/test.py | 4 +- tests/test_duration_generator.py | 110 ++++++++++++++++++ tests/test_parsing.py | 6 +- tests/test_random_generator.py | 83 +++++--------- tests/test_ts_generator.py | 17 ++- tests/test_unit.py | 190 +++++++++++++++++++++---------- 14 files changed, 519 insertions(+), 262 deletions(-) create mode 100644 src/duration_generator.py create mode 100644 src/random_generator.py create mode 100644 tests/test_duration_generator.py diff --git a/requirements.txt b/requirements.txt index 41a18c1..e1e0e43 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ numpy==1.24.4 pandas +pydantic PyYAML diff --git a/src/cluster_import.py b/src/cluster_import.py index e852d75..8bc0404 100644 --- a/src/cluster_import.py +++ b/src/cluster_import.py @@ -14,11 +14,14 @@ import numpy as np -from ts_generator import ProbilityLaw, ThermalCluster +from ts_generator import ProbabilityLaw, ThermalCluster def import_thermal_cluster(path: Path, days_per_year: int = 365) -> ThermalCluster: - law_dict = {"UNIFORM": ProbilityLaw.UNIFORM, "GEOMETRIC": ProbilityLaw.GEOMETRIC} + law_dict = { + "UNIFORM": ProbabilityLaw.UNIFORM, + "GEOMETRIC": ProbabilityLaw.GEOMETRIC, + } array = np.genfromtxt(path, delimiter=",", dtype=str) return ThermalCluster( unit_count=int(array[1][1]), @@ -28,10 +31,10 @@ def import_thermal_cluster(path: Path, days_per_year: int = 365) -> ThermalClust fo_volatility=float(array[5][1]), po_law=law_dict[array[6][1]], po_volatility=float(array[7][1]), - fo_duration=array[8][1 : days_per_year + 1].astype(float), + fo_duration=array[8][1 : days_per_year + 1].astype(int), fo_rate=array[9][1 : days_per_year + 1].astype(float), - po_duration=array[10][1 : days_per_year + 1].astype(float), + po_duration=array[10][1 : days_per_year + 1].astype(int), po_rate=array[11][1 : days_per_year + 1].astype(float), - npo_min=array[12][1 : days_per_year + 1].astype(float), - npo_max=array[13][1 : days_per_year + 1].astype(float), + npo_min=array[12][1 : days_per_year + 1].astype(int), + npo_max=array[13][1 : days_per_year + 1].astype(int), ) diff --git a/src/cluster_resolve.py b/src/cluster_resolve.py index 7ddb227..65b88c4 100644 --- a/src/cluster_resolve.py +++ b/src/cluster_resolve.py @@ -13,7 +13,7 @@ import pandas as pd from cluster_parsing import InputCluster -from ts_generator import ProbilityLaw, ThermalCluster +from ts_generator import ProbabilityLaw, ThermalCluster def resolve_thermal_cluster( @@ -21,7 +21,10 @@ def resolve_thermal_cluster( parameters_ts: pd.core.frame.DataFrame, modulation: pd.core.frame.DataFrame, ) -> ThermalCluster: - law_dict = {"UNIFORM": ProbilityLaw.UNIFORM, "GEOMETRIC": ProbilityLaw.GEOMETRIC} + law_dict = { + "UNIFORM": ProbabilityLaw.UNIFORM, + "GEOMETRIC": ProbabilityLaw.GEOMETRIC, + } return ThermalCluster( unit_count=parsed_yaml.unit_count, nominal_power=parsed_yaml.nominal_power, diff --git a/src/duration_generator.py b/src/duration_generator.py new file mode 100644 index 0000000..f330ec8 --- /dev/null +++ b/src/duration_generator.py @@ -0,0 +1,115 @@ +# Copyright (c) 2024, RTE (https://www.rte-france.com) +# +# See AUTHORS.txt +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# SPDX-License-Identifier: MPL-2.0 +# +# This file is part of the Antares project. +from abc import ABC, abstractmethod +from enum import Enum +from math import log, sqrt +from typing import List + +import numpy as np + +from random_generator import RNG + + +class ProbabilityLaw(Enum): + UNIFORM = "UNIFORM" + GEOMETRIC = "GEOMETRIC" + + +class DurationGenerator(ABC): + """ + Logic to generate unavailability duration. + """ + + @abstractmethod + def generate_duration(self, day: int) -> int: + ... + + +class GeneratorWrapper(DurationGenerator): + """ + Wraps another generator to skip unnecessary random number generations. + Used to keep backward compat with cpp implementation. + """ + + def __init__( + self, delegate: DurationGenerator, volatility: float, expecs: List[int] + ) -> None: + self.volatility = volatility + self.expectations = expecs + self.delegate = delegate + + def generate_duration(self, day: int) -> int: + """ + generation of random outage duration + """ + expectation = self.expectations[day] + # Logic copied from cpp implementation for results preservation + if self.volatility == 0 or expectation == 1: + return expectation + return self.delegate.generate_duration(day) + + +class UniformDurationGenerator(DurationGenerator): + def __init__(self, rng: RNG, volatility: float, expecs: List[int]) -> None: + self.rng = rng + self.a = np.empty(len(expecs), dtype=float) + self.b = np.empty(len(expecs), dtype=float) + for day, expec in enumerate(expecs): + xtemp = volatility * (expec - 1) + self.a[day] = expec - xtemp + self.b[day] = 2 * xtemp + 1 + + def generate_duration(self, day: int) -> int: + """ + generation of random outage duration + """ + rnd_nb = self.rng.next() + return int(self.a[day] + rnd_nb * self.b[day]) + + +class GeometricDurationGenerator(DurationGenerator): + def __init__(self, rng: RNG, volatility: float, expecs: List[int]) -> None: + self.rng = rng + self.a = np.empty(len(expecs), dtype=float) + self.b = np.empty(len(expecs), dtype=float) + for day, expec in enumerate(expecs): + xtemp = volatility * volatility * expec * (expec - 1) + if xtemp != 0: + ytemp = (sqrt(4 * xtemp + 1) - 1) / (2 * xtemp) + self.a[day] = expec - 1 / ytemp + self.b[day] = 1 / log(1 - ytemp) + else: + self.a[day] = expec - 1 + self.b[day] = 0 + + def generate_duration(self, day: int) -> int: + """ + generation of random outage duration + """ + rnd_nb = self.rng.next() + return min(int(1 + self.a[day] + self.b[day] * log(rnd_nb)), 1999) + + +def make_duration_generator( + rng: RNG, law: ProbabilityLaw, volatility: float, expectations: List[int] +) -> DurationGenerator: + """ + return a DurationGenerator for the given law + """ + base_rng: DurationGenerator + if law == ProbabilityLaw.UNIFORM: + base_rng = UniformDurationGenerator(rng, volatility, expectations) + elif law == ProbabilityLaw.GEOMETRIC: + base_rng = GeometricDurationGenerator(rng, volatility, expectations) + else: + raise ValueError(f"Unknown law type: {law}") + return GeneratorWrapper(base_rng, volatility, expectations) diff --git a/src/mersenne_twister.py b/src/mersenne_twister.py index 1976703..879ecfb 100644 --- a/src/mersenne_twister.py +++ b/src/mersenne_twister.py @@ -27,10 +27,12 @@ # LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING # NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - +import dataclasses +from dataclasses import field from typing import List, Tuple +@dataclasses.dataclass class MersenneTwister: periodN: int = 624 periodM: int = 397 @@ -40,11 +42,10 @@ class MersenneTwister: MAG: Tuple[int, int] = (0, MATRIX_A) - mt: List[int] = [0] * periodN + mt: List[int] = field(default_factory=lambda: [0] * 624) mti: int = 0 - @classmethod def seed(self, seed: int) -> None: self.mt[0] = seed & 0xFFFFFFFF for i in range(1, self.periodN): @@ -53,7 +54,6 @@ def seed(self, seed: int) -> None: self.mti = self.periodN - @classmethod def next(self) -> float: if self.mti == self.periodN: for j in range(self.periodN - self.periodM): @@ -87,6 +87,5 @@ def next(self) -> float: return y / 4294967295 - @classmethod def reset(self) -> None: self.seed(5489) diff --git a/src/random_generator.py b/src/random_generator.py new file mode 100644 index 0000000..dc78072 --- /dev/null +++ b/src/random_generator.py @@ -0,0 +1,48 @@ +# Copyright (c) 2024, RTE (https://www.rte-france.com) +# +# See AUTHORS.txt +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# SPDX-License-Identifier: MPL-2.0 +# +# This file is part of the Antares project. +import random +from abc import ABC, abstractmethod +from typing import Optional + +from mersenne_twister import MersenneTwister + + +class RNG(ABC): + """ + Random number generator interface + """ + + @abstractmethod + def next(self) -> float: + ... + + +class PythonRNG(ABC): + """ + Native python RNG. + """ + + def next(self) -> float: + return random.random() + + +class MersenneTwisterRNG(RNG): + """ + Our own RNG based on Mersenne-Twister algorithm. + """ + + def __init__(self, seed: int = 5489): + self._rng = MersenneTwister() + self._rng.seed(seed) + + def next(self) -> float: + return self._rng.next() diff --git a/src/ts_generator.py b/src/ts_generator.py index 6c357fe..59eb93b 100644 --- a/src/ts_generator.py +++ b/src/ts_generator.py @@ -10,86 +10,19 @@ # # This file is part of the Antares project. -from abc import ABC, abstractmethod from dataclasses import dataclass -from enum import Enum -from math import log, sqrt -from random import random -from typing import Any, List, Optional, Tuple +from typing import List import numpy as np -from numpy.typing import ArrayLike +import pandas as pd -from mersenne_twister import MersenneTwister +from duration_generator import ProbabilityLaw, make_duration_generator +from random_generator import RNG, MersenneTwisterRNG # probabilities above FAILURE_RATE_EQ_1 are considered certain (equal to 1) FAILURE_RATE_EQ_1 = 0.999 -class PythonGenerator: - @classmethod - def next(self) -> float: - return random() - - -rndgenerator = MersenneTwister -rndgenerator.seed(1) - - -class ProbilityLaw(Enum): - UNIFORM = "UNIFORM" - GEOMETRIC = "GEOMETRIC" - - -class DurationGenerator(ABC): - @abstractmethod - def __init__(self, volatility: float, expecs: List[int]) -> None: - ... - - @abstractmethod - def generate_duration(self, day: int) -> int: - ... - - -class UniformDurationGenerator(DurationGenerator): - def __init__(self, volatility: float, expecs: List[int]) -> None: - self.a = np.empty(len(expecs), dtype=float) - self.b = np.empty(len(expecs), dtype=float) - for day, expec in enumerate(expecs): - xtemp = volatility * (expec - 1) - self.a[day] = expec - xtemp - self.b[day] = 2 * xtemp + 1 - - def generate_duration(self, day: int) -> int: - """ - generation of random outage duration - """ - rnd_nb = rndgenerator.next() - return int(self.a[day] + rnd_nb * self.b[day]) - - -class GeometricDurationGenerator(DurationGenerator): - def __init__(self, volatility: float, expecs: List[int]) -> None: - self.a = np.empty(len(expecs), dtype=float) - self.b = np.empty(len(expecs), dtype=float) - for day, expec in enumerate(expecs): - xtemp = volatility * volatility * expec * (expec - 1) - if xtemp != 0: - ytemp = (sqrt(4 * xtemp + 1) - 1) / (2 * xtemp) - self.a[day] = expec - 1 / ytemp - self.b[day] = 1 / log(1 - ytemp) - else: - self.a[day] = expec - 1 - self.b[day] = 0 - - def generate_duration(self, day: int) -> int: - """ - generation of random outage duration - """ - rnd_nb = rndgenerator.next() - return min(int(1 + self.a[day] + self.b[day] * log(rnd_nb)), 1999) - - @dataclass class ThermalCluster: # available units of the cluster @@ -97,6 +30,7 @@ class ThermalCluster: # nominal power nominal_power: float # modulation of the nominal power for a certain hour in the day (between 0 and 1) + # TODO: check that it should be 24 or 8760 ? modulation: List[float] ### maybe group nominal_power and modulation in one vaiable # forced and planed outage parameters @@ -110,30 +44,33 @@ class ThermalCluster: # forced and planed outage probability law and volatility # volatility characterizes the distance from the expect at which the value drawn can be - fo_law: ProbilityLaw + fo_law: ProbabilityLaw fo_volatility: float - po_law: ProbilityLaw + po_law: ProbabilityLaw po_volatility: float class OutputTimeseries: - def __init__(self, ts_nb: int, days_per_year: int) -> None: + def __init__(self, ts_count: int, days: int) -> None: + self.available_units = np.zeros(shape=(ts_count, days), dtype=int) # available power each hours - self.available_power = np.empty((ts_nb, 24 * days_per_year), dtype=float) + self.available_power = np.zeros((ts_count, 24 * days), dtype=float) # number of pure planed, pure forced and mixed outage each day - self.nb_ppo = np.empty((ts_nb, days_per_year), dtype=int) - self.nb_pfo = np.empty((ts_nb, days_per_year), dtype=int) - self.nb_mxo = np.empty((ts_nb, days_per_year), dtype=int) + self.planned_outages = np.zeros((ts_count, days), dtype=int) + self.forced_outages = np.zeros((ts_count, days), dtype=int) + self.mixed_outages = np.zeros((ts_count, days), dtype=int) # number of pure planed and pure forced outage duration each day # (mixed outage duration = pod + fod) - self.pod = np.empty((ts_nb, days_per_year), dtype=int) - self.fod = np.empty((ts_nb, days_per_year), dtype=int) + self.planned_outage_durations = np.zeros((ts_count, days), dtype=int) + self.forced_outage_durations = np.zeros((ts_count, days), dtype=int) class ThermalDataGenerator: - def __init__(self, days_per_year: int = 365) -> None: - self.days_per_year = days_per_year + def __init__(self, rng: RNG = MersenneTwisterRNG(), days: int = 365) -> None: + self.rng = rng + self.days = days + # TODO: Remove this log size limits, seems useless and error prone if very large durations self.log_size = 4000 # >= 5 * (max(df) + max(dp)) # the number of starting (if positive)/ stopping (if negative) units (due to FO and PO) at a given time self.LOG = [0] * self.log_size @@ -141,23 +78,12 @@ def __init__(self, days_per_year: int = 365) -> None: self.LOGP = [0] * self.log_size # pogramed and forced outage rate - self.lf = np.empty(days_per_year, dtype=float) - self.lp = np.empty(days_per_year, dtype=float) + self.lf = np.zeros(days, dtype=float) + self.lp = np.zeros(days, dtype=float) ## ??? - self.ff = np.empty(days_per_year, dtype=float) # ff = lf / (1 - lf) - self.pp = np.empty(days_per_year, dtype=float) # pp = lp / (1 - lp) - - def get_duration_generator( - self, law: ProbilityLaw, volatility: float, expecs: List[int] - ) -> DurationGenerator: - """ - return a DurationGenerator for the given law - """ - if law == ProbilityLaw.UNIFORM: - return UniformDurationGenerator(volatility, expecs) - elif law == ProbilityLaw.GEOMETRIC: - return GeometricDurationGenerator(volatility, expecs) + self.ff = np.zeros(days, dtype=float) # ff = lf / (1 - lf) + self.pp = np.zeros(days, dtype=float) # pp = lp / (1 - lp) def generate_time_series( self, @@ -173,7 +99,7 @@ def generate_time_series( self.FPOW: List[List[float]] = [] self.PPOW: List[List[float]] = [] - for day in range(self.days_per_year): + for day in range(self.days): # lf and lp represent the forced and programed failure rate # failure rate means the probability to enter in outage each day # its value is given by: OR / [OR + OD * (1 - OR)] @@ -213,11 +139,11 @@ def generate_time_series( self.FPOW[-1].append(pow(a, k)) self.PPOW[-1].append(pow(b, k)) - self.fod_generator = self.get_duration_generator( - cluster.fo_law, cluster.fo_volatility, cluster.fo_duration + fod_generator = make_duration_generator( + self.rng, cluster.fo_law, cluster.fo_volatility, cluster.fo_duration ) - self.pod_generator = self.get_duration_generator( - cluster.po_law, cluster.po_volatility, cluster.po_duration + pod_generator = make_duration_generator( + self.rng, cluster.po_law, cluster.po_volatility, cluster.po_duration ) # --- calculation --- @@ -225,7 +151,7 @@ def generate_time_series( # as a consequence, N + 2 time series will be computed # output that will be returned - output = OutputTimeseries(number_of_timeseries, self.days_per_year) + output = OutputTimeseries(number_of_timeseries, self.days) # mixed, pure planned and pure force outage MXO = 0 @@ -247,7 +173,7 @@ def generate_time_series( # hour in the year hour = 0 - for day in range(self.days_per_year): + for day in range(self.days): # = return of units wich were in outage = cur_nb_PO -= self.LOGP[now] self.LOGP[ @@ -256,13 +182,11 @@ def generate_time_series( cur_nb_AU += self.LOG[now] self.LOG[now] = 0 - # = determinating units that go on outage = - # FO and PO canditate FOC = 0 POC = 0 if self.lf[day] > 0 and self.lf[day] <= FAILURE_RATE_EQ_1: - A = rndgenerator.next() + A = self.rng.next() last = self.FPOW[day][cur_nb_AU] if A > last: cumul = last @@ -272,7 +196,9 @@ def generate_time_series( FOC = d if A <= cumul: break - elif self.lf[day] > FAILURE_RATE_EQ_1: + elif ( + self.lf[day] > FAILURE_RATE_EQ_1 + ): # TODO: not same comparison as cpp ? FOC = cur_nb_AU else: # self.lf[day] == 0 FOC = 0 @@ -285,7 +211,7 @@ def generate_time_series( elif stock > cur_nb_AU: AUN_app = 0 - A = rndgenerator.next() + A = self.rng.next() last = self.PPOW[day][cur_nb_AU] if A > last: cumul = last @@ -359,9 +285,9 @@ def generate_time_series( true_FOD = 0 if PPO != 0 or MXO != 0: - true_POD = self.pod_generator.generate_duration(day) + true_POD = pod_generator.generate_duration(day) if PFO != 0 or MXO != 0: - true_FOD = self.fod_generator.generate_duration(day) + true_FOD = fod_generator.generate_duration(day) if PPO != 0: fut = (now + true_POD) % self.log_size @@ -377,17 +303,18 @@ def generate_time_series( # = storing output in output arrays = if ts_index >= 0: # drop the 2 first generated timeseries - output.nb_ppo[ts_index][day] = PPO - output.nb_pfo[ts_index][day] = PFO - output.nb_mxo[ts_index][day] = MXO - output.pod[ts_index][day] = true_POD - output.fod[ts_index][day] = true_FOD - for h in range(24): - output.available_power[ts_index][hour] = ( - cur_nb_AU * cluster.nominal_power * cluster.modulation[h] - ) - hour += 1 + output.planned_outages[ts_index][day] = PPO + output.forced_outages[ts_index][day] = PFO + output.mixed_outages[ts_index][day] = MXO + output.planned_outage_durations[ts_index][day] = true_POD + output.forced_outage_durations[ts_index][day] = true_FOD + output.available_units[ts_index][day] = cur_nb_AU now = (now + 1) % self.log_size + output.available_power = ( + np.repeat(output.available_units, 24, axis=1) + * cluster.nominal_power + * np.tile(cluster.modulation, self.days) + ) return output diff --git a/tests/conftest.py b/tests/conftest.py index d9d3ece..3dd5726 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -14,6 +14,8 @@ import pytest +from random_generator import RNG, MersenneTwisterRNG + @pytest.fixture def data_directory() -> Path: @@ -23,3 +25,8 @@ def data_directory() -> Path: @pytest.fixture def output_directory() -> Path: return Path(__file__).parent / "output" + + +@pytest.fixture +def rng() -> RNG: + return MersenneTwisterRNG() diff --git a/tests/test.py b/tests/test.py index 970d291..4dc125b 100644 --- a/tests/test.py +++ b/tests/test.py @@ -33,8 +33,8 @@ def test_cluster(cluster, output_directory): tot_po = 0 tot_fo = 0 for i in range(365 * ts_nb): - tot_po += results.nb_ppo[i // 365][i % 365] * 2 - tot_fo += results.nb_pfo[i // 365][i % 365] * 8 + tot_po += results.planned_outages[i // 365][i % 365] * 2 + tot_fo += results.forced_outages[i // 365][i % 365] * 8 true_por = tot_po / (365 * ts_nb) true_for = tot_fo / (365 * ts_nb) diff --git a/tests/test_duration_generator.py b/tests/test_duration_generator.py new file mode 100644 index 0000000..e0bd3d4 --- /dev/null +++ b/tests/test_duration_generator.py @@ -0,0 +1,110 @@ +# Copyright (c) 2024, RTE (https://www.rte-france.com) +# +# See AUTHORS.txt +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# SPDX-License-Identifier: MPL-2.0 +# +# This file is part of the Antares project. +from unittest import expectedFailure + +import matplotlib.pyplot as plt +import pytest + +from duration_generator import ( + GeometricDurationGenerator, + ProbabilityLaw, + UniformDurationGenerator, + make_duration_generator, +) +from random_generator import MersenneTwisterRNG + + +@pytest.mark.parametrize( + "volatility,expectations,expected", + [ + (0, [4, 2, 3, 1], [4, 2, 3, 1]), # Expecting output = input + (1, [4, 2, 3, 1], [6, 1, 5, 1]), + (2, [4, 2, 3, 1], [8, 0, 7, 1]), # Expecting twice the deviation as above + ], +) +def test_uniform_law_generator(rng, volatility, expectations, expected): + generator = UniformDurationGenerator(rng, volatility, expectations) + assert [generator.generate_duration(d) for d in range(len(expected))] == expected + + +# TODO: results are suprisingly low, check this. Negative results can be obtained... +@pytest.mark.parametrize( + "volatility,expectations,expected", + [ + (0, [4, 2, 3, 1], [4, 2, 3, 1]), # Expecting output = input + (1, [4, 2, 3, 1], [1, 3, 1, 1]), + (2, [4, 2, 3, 1], [-1, 5, 0, 1]), # Expecting twice the deviation as above + ], +) +def test_geometric_law_generator(rng, volatility, expectations, expected): + generator = GeometricDurationGenerator(rng, volatility, expectations) + assert [generator.generate_duration(d) for d in range(len(expected))] == expected + + +def test_legacy_generator_skips_rng_when_zero_vol(): + rng = MersenneTwisterRNG() + generator = make_duration_generator( + rng, ProbabilityLaw.UNIFORM, volatility=0, expectations=[10, 10] + ) + generator.generate_duration(0) + # Check we still have a random number identical to the first one that should be generated + assert rng.next() == MersenneTwisterRNG().next() + + +def test_geometric_law(rng, output_directory): + volatility = 1 + generator = GeometricDurationGenerator(rng, volatility, [10]) + + expec = 0 + nb_values = 45 + values = [0] * nb_values + N = 1000000 + N_inv = 1 / N + for _ in range(N): + value = generator.generate_duration(0) + assert value >= 1 + expec += value + + if value < nb_values: + values[value] += N_inv + + expec /= N + assert expec == pytest.approx(10, abs=0.1) + + plt.plot(range(nb_values), values) + plt.savefig(output_directory / "geometric_law_distrib.png") + plt.clf() + + +def test_uniform_law(rng, output_directory): + volatility = 1 + generator = UniformDurationGenerator(rng, volatility, [10]) + + expec = 0 + nb_values = 45 + values = [0] * nb_values + N = 1000000 + N_inv = 1 / N + for _ in range(N): + value = generator.generate_duration(0) + assert value >= 1 + expec += value + + if value < nb_values: + values[value] += N_inv + + expec /= N + assert expec == pytest.approx(10, abs=0.1) + + plt.plot(range(nb_values), values) + plt.savefig(output_directory / "uniform_law_distrib.png") + plt.clf() diff --git a/tests/test_parsing.py b/tests/test_parsing.py index dd9b413..c89e9fd 100644 --- a/tests/test_parsing.py +++ b/tests/test_parsing.py @@ -14,7 +14,7 @@ from cluster_parsing import parse_cluster_ts, parse_yaml_cluster from cluster_resolve import resolve_thermal_cluster -from ts_generator import ProbilityLaw, ThermalCluster, ThermalDataGenerator +from ts_generator import ProbabilityLaw, ThermalCluster, ThermalDataGenerator NB_OF_DAY = 10 @@ -25,9 +25,9 @@ def cluster() -> ThermalCluster: unit_count=1, nominal_power=500, modulation=[1 for i in range(24)], - fo_law=ProbilityLaw.UNIFORM, + fo_law=ProbabilityLaw.UNIFORM, fo_volatility=0, - po_law=ProbilityLaw.UNIFORM, + po_law=ProbabilityLaw.UNIFORM, po_volatility=0, fo_duration=[2 for i in range(NB_OF_DAY)], fo_rate=[0.2 for i in range(NB_OF_DAY)], diff --git a/tests/test_random_generator.py b/tests/test_random_generator.py index 6faaf59..928deca 100644 --- a/tests/test_random_generator.py +++ b/tests/test_random_generator.py @@ -9,59 +9,30 @@ # SPDX-License-Identifier: MPL-2.0 # # This file is part of the Antares project. - -import matplotlib.pyplot as plt -import pytest - -from cluster_import import import_thermal_cluster -from ts_generator import GeometricDurationGenerator, UniformDurationGenerator - - -def test_geometric_law(output_directory): - volatility = 1 - generator = GeometricDurationGenerator(volatility, [10]) - - expec = 0 - nb_values = 45 - values = [0] * nb_values - N = 1000000 - N_inv = 1 / N - for _ in range(N): - value = generator.generate_duration(0) - assert value >= 1 - expec += value - - if value < nb_values: - values[value] += N_inv - - expec /= N - assert expec == pytest.approx(10, abs=0.1) - - plt.plot(range(nb_values), values) - plt.savefig(output_directory / "geometric_law_distrib.png") - plt.clf() - - -def test_uniform_law(output_directory): - volatility = 1 - generator = UniformDurationGenerator(volatility, [10]) - - expec = 0 - nb_values = 45 - values = [0] * nb_values - N = 1000000 - N_inv = 1 / N - for _ in range(N): - value = generator.generate_duration(0) - assert value >= 1 - expec += value - - if value < nb_values: - values[value] += N_inv - - expec /= N - assert expec == pytest.approx(10, abs=0.1) - - plt.plot(range(nb_values), values) - plt.savefig(output_directory / "uniform_law_distrib.png") - plt.clf() +import random + +from random_generator import MersenneTwisterRNG, PythonRNG + + +def test_python_rng(): + rng = PythonRNG() + random.seed(10) + expected = [random.random() for _ in range(10)] + random.seed(10) + assert [rng.next() for _ in range(10)] == expected + + +def test_custom_rng(): + rng = MersenneTwisterRNG() + assert [rng.next() for _ in range(10)] == [ + 0.8147236920927473, + 0.13547700413863104, + 0.9057919343248456, + 0.835008589978099, + 0.12698681189841285, + 0.968867771320247, + 0.9133758558690026, + 0.22103404282150652, + 0.6323592501302154, + 0.30816705043152137, + ] diff --git a/tests/test_ts_generator.py b/tests/test_ts_generator.py index f36ea69..bb399c4 100644 --- a/tests/test_ts_generator.py +++ b/tests/test_ts_generator.py @@ -11,7 +11,6 @@ # This file is part of the Antares project. import csv -from pathlib import Path import pytest @@ -43,8 +42,8 @@ def test_one_unit_cluster(cluster_1, output_directory): tot_po = 0 tot_fo = 0 for i in range(365 * ts_nb): - tot_po += results.nb_ppo[i // 365][i % 365] * 2 - tot_fo += results.nb_pfo[i // 365][i % 365] * 8 + tot_po += results.planned_outages[i // 365][i % 365] * 2 + tot_fo += results.forced_outages[i // 365][i % 365] * 8 true_por = tot_po / (365 * ts_nb) true_for = tot_fo / (365 * ts_nb) @@ -74,8 +73,8 @@ def test_hundred_unit_cluster(cluster_100, output_directory): tot_po = 0 tot_fo = 0 for i in range(365 * ts_nb): - tot_po += results.nb_ppo[i // 365][i % 365] * 2 - tot_fo += results.nb_pfo[i // 365][i % 365] * 8 + tot_po += results.planned_outages[i // 365][i % 365] * 2 + tot_fo += results.forced_outages[i // 365][i % 365] * 8 true_por = tot_po / (365 * ts_nb) true_for = tot_fo / (365 * ts_nb) @@ -84,8 +83,8 @@ def test_hundred_unit_cluster(cluster_100, output_directory): cursor = [0] * 10 tot_simult_po = 0 for i in range(365 * ts_nb): - po = results.nb_ppo[i // 365][i % 365] - mo = results.nb_mxo[i // 365][i % 365] + po = results.planned_outages[i // 365][i % 365] + mo = results.mixed_outages[i // 365][i % 365] tot_simult_po += po tot_simult_po += mo @@ -132,8 +131,8 @@ def test_max_po(cluster_high_por, output_directory): cursor = [0] * 10 tot_simult_po = 0 for i in range(365 * ts_nb): - po = results.nb_ppo[i // 365][i % 365] - mo = results.nb_mxo[i // 365][i % 365] + po = results.planned_outages[i // 365][i % 365] + mo = results.mixed_outages[i // 365][i % 365] tot_simult_po += po tot_simult_po += mo diff --git a/tests/test_unit.py b/tests/test_unit.py index c26c5e6..863874b 100644 --- a/tests/test_unit.py +++ b/tests/test_unit.py @@ -9,71 +9,145 @@ # SPDX-License-Identifier: MPL-2.0 # # This file is part of the Antares project. +import cProfile +import random +from pstats import SortKey -import pytest +import numpy as np +import numpy.testing as npt -from ts_generator import ProbilityLaw, ThermalCluster, ThermalDataGenerator +from ts_generator import ProbabilityLaw, ThermalCluster, ThermalDataGenerator -NB_OF_DAY = 10 - -@pytest.fixture -def cluster() -> ThermalCluster: - return ThermalCluster( - unit_count=1, - nominal_power=500, +def test_forced_outages(rng): + days = 365 + cluster = ThermalCluster( + unit_count=10, + nominal_power=100, modulation=[1 for i in range(24)], - fo_law=ProbilityLaw.UNIFORM, + fo_law=ProbabilityLaw.UNIFORM, fo_volatility=0, - po_law=ProbilityLaw.UNIFORM, + po_law=ProbabilityLaw.UNIFORM, po_volatility=0, - fo_duration=[2 for i in range(NB_OF_DAY)], - fo_rate=[0.2 for i in range(NB_OF_DAY)], - po_duration=[1 for i in range(NB_OF_DAY)], - po_rate=[0.1 for i in range(NB_OF_DAY)], - npo_min=[0 for i in range(NB_OF_DAY)], - npo_max=[1 for i in range(NB_OF_DAY)], + fo_duration=[10 for i in range(days)], + fo_rate=[0.2 for i in range(days)], + po_duration=[10 for i in range(days)], + po_rate=[0 for i in range(days)], + npo_min=[0 for i in range(days)], + npo_max=[10 for i in range(days)], + ) + cluster.modulation[12] = 0.5 + + generator = ThermalDataGenerator(rng=rng, days=days) + results = generator.generate_time_series(cluster, 1) + # 2 forced outages occur on day 5, with duration 10 + npt.assert_equal(results.forced_outages[0][:6], [0, 0, 0, 0, 2, 0]) + npt.assert_equal(results.forced_outage_durations[0][:6], [0, 0, 0, 0, 10, 0]) + # No planned outage + npt.assert_equal(results.planned_outages[0], np.zeros(365)) + npt.assert_equal(results.planned_outage_durations[0], np.zeros(365)) + + npt.assert_equal(results.available_units[0][:5], [9, 9, 9, 9, 8]) + # Check available power consistency with available units and modulation + assert results.available_power[0][0] == 900 + assert results.available_power[0][12] == 450 # Modulation is 0.5 for hour 12 + assert results.available_power[0][4 * 24] == 800 + + +def test_planned_outages(rng): + days = 365 + cluster = ThermalCluster( + unit_count=10, + nominal_power=100, + modulation=[1 for i in range(24)], + fo_law=ProbabilityLaw.UNIFORM, + fo_volatility=0, + po_law=ProbabilityLaw.UNIFORM, + po_volatility=0, + fo_duration=[10 for i in range(days)], + fo_rate=[0 for i in range(days)], + po_duration=[10 for i in range(days)], + po_rate=[0.2 for i in range(days)], + npo_min=[0 for i in range(days)], + npo_max=[10 for i in range(days)], + ) + cluster.modulation[12] = 0.5 + + generator = ThermalDataGenerator(rng=rng, days=days) + results = generator.generate_time_series(cluster, 1) + # 0 forced outage + npt.assert_equal(results.forced_outages[0], np.zeros(365)) + npt.assert_equal(results.forced_outage_durations[0], np.zeros(365)) + # No planned outage + npt.assert_equal(results.planned_outages[0][:6], [0, 0, 0, 0, 2, 0]) + npt.assert_equal(results.available_units[0][:5], [9, 9, 9, 9, 8]) + # Check available power consistency with available units and modulation + assert results.available_power[0][0] == 900 + assert results.available_power[0][12] == 450 # Modulation is 0.5 for hour 12 + assert results.available_power[0][4 * 24] == 800 + + +def test_planned_outages_limitation(rng): + days = 365 + # Maximum 1 planned outage at a time. + cluster = ThermalCluster( + unit_count=10, + nominal_power=100, + modulation=[1 for i in range(24)], + fo_law=ProbabilityLaw.UNIFORM, + fo_volatility=0, + po_law=ProbabilityLaw.UNIFORM, + po_volatility=0, + fo_duration=[10 for i in range(days)], + fo_rate=[0 for i in range(days)], + po_duration=[2 for i in range(days)], + po_rate=[0.2 for i in range(days)], + npo_min=[0 for i in range(days)], + npo_max=[1 for i in range(days)], ) + generator = ThermalDataGenerator(rng=rng, days=days) + results = generator.generate_time_series(cluster, 1) + # No forced outage + npt.assert_equal(results.forced_outages[0], np.zeros(365)) + npt.assert_equal(results.forced_outage_durations[0], np.zeros(365)) + # Maxmimum one planned outage at a time + npt.assert_equal(results.planned_outages[0][:6], [0, 1, 0, 1, 0, 1]) + npt.assert_equal(results.planned_outage_durations[0][:6], [0, 2, 0, 2, 0, 2]) + npt.assert_equal(results.available_units[0][:5], [9, 9, 9, 9, 9]) + # Check available power consistency with available units and modulation + assert results.available_power[0][0] == 900 + assert results.available_power[0][4 * 24] == 900 + + +def test_planned_outages_min_limitation(rng): + days = 365 + # Minimum 2 planned outages at a time + cluster = ThermalCluster( + unit_count=10, + nominal_power=100, + modulation=[1 for i in range(24)], + fo_law=ProbabilityLaw.UNIFORM, + fo_volatility=0, + po_law=ProbabilityLaw.UNIFORM, + po_volatility=0, + fo_duration=[10 for i in range(days)], + fo_rate=[0 for i in range(days)], + po_duration=[10 for i in range(days)], + po_rate=[0.2 for i in range(days)], + npo_min=[2 for i in range(days)], + npo_max=[5 for i in range(days)], + ) -def test_ts_value(cluster): - ts_nb = 4 - - generator = ThermalDataGenerator(days_per_year=NB_OF_DAY) - results = generator.generate_time_series(cluster, ts_nb) - - assert results.available_power.shape == (ts_nb, NB_OF_DAY * 24) - assert results.nb_ppo.shape == (ts_nb, NB_OF_DAY) - assert results.nb_pfo.shape == (ts_nb, NB_OF_DAY) - assert results.nb_mxo.shape == (ts_nb, NB_OF_DAY) - assert results.pod.shape == (ts_nb, NB_OF_DAY) - assert results.fod.shape == (ts_nb, NB_OF_DAY) - - for l in range(ts_nb): - for c in range(NB_OF_DAY * 24): - assert results.available_power[l][c] % 500 == 0 - - -def test_ts_value_with_modulation(cluster): - modulation = [(i % 10) * 0.1 for i in range(24)] - - cluster.modulation = modulation - - ts_nb = 4 - - generator = ThermalDataGenerator(days_per_year=NB_OF_DAY) - results = generator.generate_time_series(cluster, ts_nb) - - assert results.available_power.shape == (ts_nb, NB_OF_DAY * 24) - assert results.nb_ppo.shape == (ts_nb, NB_OF_DAY) - assert results.nb_pfo.shape == (ts_nb, NB_OF_DAY) - assert results.nb_mxo.shape == (ts_nb, NB_OF_DAY) - assert results.pod.shape == (ts_nb, NB_OF_DAY) - assert results.fod.shape == (ts_nb, NB_OF_DAY) - - for l in range(ts_nb): - for d in range(NB_OF_DAY): - if modulation[d] != 0: - assert results.available_power[l][d * 24] % (500 * modulation[d]) == 0 - else: - assert results.available_power[l][d * 24] == 0 + generator = ThermalDataGenerator(rng=rng, days=days) + results = generator.generate_time_series(cluster, 1) + # No forced outage + npt.assert_equal(results.forced_outages[0], np.zeros(365)) + npt.assert_equal(results.forced_outage_durations[0], np.zeros(365)) + # Maxmimum one planned outage at a time + npt.assert_equal(results.planned_outages[0][:6], [0, 0, 1, 0, 0, 1]) + npt.assert_equal(results.planned_outage_durations[0][:6], [0, 0, 10, 0, 0, 10]) + npt.assert_equal(results.available_units[0][:5], [8, 8, 8, 8, 8]) + # Check available power consistency with available units and modulation + assert results.available_power[0][0] == 800 + assert results.available_power[0][4 * 24] == 800