Hydro heuristic

src/andromede/hydro_heuristic/data.py

@@ -0,0 +1,175 @@
+# 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 dataclasses import dataclass
+from pathlib import Path
+from typing import List, Optional
+
+import numpy as np
+
+
+@dataclass
+class ReservoirParameters:
+    capacity: float
+    initial_level: float
+    folder_name: str
+    scenario: int
+
+
+@dataclass
+class HydroHeuristicParameters:
+    inter_breakdown: int = 1
+    total_target: Optional[float] = None
+
+
+@dataclass
+class DataAggregatorParameters:
+    hours_aggregated_time_steps: List[int]
+    timesteps: List[int]
+
+
+def get_number_of_days_in_month(month: int) -> int:
+    number_day_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31][month]
+    return number_day_month
+
+
+class RawHydroData:
+    def __init__(self, folder_name: str, scenario: int) -> None:
+        self.folder_name = folder_name
+        self.scenario = scenario
+
+        self.name_file = {
+            "demand": "load",
+            "inflow": "mod",
+            "lower_rule_curve": "reservoir",
+            "upper_rule_curve": "reservoir",
+            "max_generating": "maxpower",
+        }
+
+        self.column = {
+            "demand": self.scenario,
+            "inflow": self.scenario,
+            "lower_rule_curve": 0,
+            "upper_rule_curve": 2,
+            "max_generating": 0,
+        }
+
+    def read_data(self, name: str) -> list[float]:
+        hours_input = 1 if name == "demand" else 24
+
+        data = np.loadtxt(
+            Path(__file__).parent
+            / (
+                "../../../tests/functional/data/"
+                + self.folder_name
+                + "/"
+                + self.name_file[name]
+                + ".txt"
+            )
+        )
+        if len(data.shape) >= 2:
+            data = data[:, self.column[name]]
+        data = np.repeat(data, hours_input)
+        if self.name_file[name] == "mod":
+            data = data / hours_input
+
+        return list(data)
+
+
+class HydroHeuristicData:
+    def __init__(
+        self,
+        data_aggregator_parameters: DataAggregatorParameters,
+        reservoir_data: ReservoirParameters,
+    ):
+        self.reservoir_data = reservoir_data
+
+        data_aggregator = DataAggregator(
+            data_aggregator_parameters.hours_aggregated_time_steps,
+            data_aggregator_parameters.timesteps,
+        )
+
+        raw_data_reader = RawHydroData(
+            reservoir_data.folder_name, reservoir_data.scenario
+        )
+
+        self.demand = data_aggregator.aggregate_data(
+            operator="sum",
+            data=raw_data_reader.read_data("demand"),
+        )
+        self.inflow = data_aggregator.aggregate_data(
+            operator="sum",
+            data=raw_data_reader.read_data("inflow"),
+        )
+        self.lower_rule_curve = data_aggregator.aggregate_data(
+            operator="lag_first_element",
+            data=raw_data_reader.read_data("lower_rule_curve"),
+        )
+        self.upper_rule_curve = data_aggregator.aggregate_data(
+            operator="lag_first_element",
+            data=raw_data_reader.read_data("upper_rule_curve"),
+        )
+        self.max_generating = data_aggregator.aggregate_data(
+            operator="sum",
+            data=raw_data_reader.read_data("max_generating"),
+        )
+
+    def compute_target(self, heuristic_parameters: HydroHeuristicParameters) -> None:
+        if heuristic_parameters.total_target is None:
+            total_target = sum(self.inflow)
+        else:
+            total_target = heuristic_parameters.total_target
+        target = (
+            total_target
+            * np.power(self.demand, heuristic_parameters.inter_breakdown)
+            / sum(np.power(self.demand, heuristic_parameters.inter_breakdown))
+        )
+
+        self.target = list(target)
+
+
+class DataAggregator:
+    def __init__(
+        self,
+        hours_aggregated_time_steps: List[int],
+        timesteps: List[int],
+    ) -> None:
+        self.hours_aggregated_time_steps = hours_aggregated_time_steps
+        self.timesteps = timesteps
+
+    def aggregate_data(self, operator: str, data: list[float]) -> List[float]:
+        aggregated_data: List[float] = []
+        hour = 0
+        for time_step, hours_time_step in enumerate(self.hours_aggregated_time_steps):
+            if time_step in self.timesteps:
+                if operator == "sum":
+                    aggregated_data.append(np.sum(data[hour : hour + hours_time_step]))
+                elif operator == "lag_first_element":
+                    aggregated_data.append(data[(hour + hours_time_step) % len(data)])
+            hour += hours_time_step
+        return aggregated_data
+
+
+def update_generation_target(
+    all_daily_generation: list[float], daily_generation: list[float]
+) -> list[float]:
+    all_daily_generation = all_daily_generation + daily_generation
+    return all_daily_generation
+
+
+def calculate_weekly_target(all_daily_generation: list[float]) -> list[float]:
+    weekly_target = [
+        sum([all_daily_generation[day] for day in range(7 * week, 7 * (week + 1))])
+        for week in range(len(all_daily_generation) // 7)
+    ]
+
+    return weekly_target

src/andromede/hydro_heuristic/heuristic_model.py

@@ -0,0 +1,173 @@
+# 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 typing import List
+
+from andromede.expression import ExpressionNode, literal, param, var
+from andromede.expression.indexing_structure import IndexingStructure
+from andromede.model import Model, float_parameter, float_variable, model
+from andromede.model.constraint import Constraint
+from andromede.model.parameter import Parameter, float_parameter
+from andromede.model.variable import Variable, float_variable
+
+CONSTANT = IndexingStructure(False, False)
+TIME_AND_SCENARIO_FREE = IndexingStructure(True, True)
+ANTICIPATIVE_TIME_VARYING = IndexingStructure(True, True)
+NON_ANTICIPATIVE_TIME_VARYING = IndexingStructure(True, False)
+CONSTANT_PER_SCENARIO = IndexingStructure(False, True)
+
+
+class HeuristicHydroModelBuilder:
+    def __init__(
+        self,
+        hydro_model: Model,
+        horizon: str,
+    ) -> None:
+        self.hydro_model = hydro_model
+        self.horizon = horizon
+
+        self.objective_function_cost = {
+            "gamma_d": 1,
+            "gamma_delta": 1 if horizon == "monthly" else 2,
+            "gamma_y": 100000 if horizon == "monthly" else 68,
+            "gamma_w": 0 if horizon == "monthly" else 34,
+            "gamma_v+": 100 if horizon == "monthly" else 0,
+            "gamma_v-": 100 if horizon == "monthly" else 68,
+            "gamma_o": 0 if horizon == "monthly" else 23 * 68 + 1,
+            "gamma_s": 0 if horizon == "monthly" else -1 / 32,
+        }
+
+    def get_model(
+        self,
+    ) -> Model:
+        assert "level" in self.hydro_model.variables.keys()
+        assert "generating" in self.hydro_model.variables.keys()
+        assert "overflow" in self.hydro_model.variables.keys()
+        assert "lower_rule_curve" in self.hydro_model.parameters.keys()
+        assert "upper_rule_curve" in self.hydro_model.parameters.keys()
+        HYDRO_HEURISTIC = model(
+            id="H",
+            parameters=list(self.hydro_model.parameters.values())
+            + self.get_heuristic_parameters(),
+            variables=list(self.hydro_model.variables.values())
+            + self.get_heuristic_variables(),
+            constraints=list(self.hydro_model.constraints.values())
+            + self.get_heuristic_constraints(),
+            objective_operational_contribution=self.get_heuristic_objective(),
+        )
+        return HYDRO_HEURISTIC
+
+    def get_heuristic_objective(self) -> ExpressionNode:
+        return (
+            self.objective_function_cost["gamma_d"]
+            * var("distance_between_target_and_generating")
+            + self.objective_function_cost["gamma_v+"]
+            * var("violation_upper_rule_curve")
+            + self.objective_function_cost["gamma_v-"]
+            * var("violation_lower_rule_curve")
+            + self.objective_function_cost["gamma_o"] * var("overflow")
+            + self.objective_function_cost["gamma_s"] * var("level")
+        ).sum().expec() + (
+            self.objective_function_cost["gamma_delta"]
+            * var("max_distance_between_target_and_generating")
+            + self.objective_function_cost["gamma_y"]
+            * var("max_violation_lower_rule_curve")
+            + self.objective_function_cost["gamma_w"] * var("gap_to_target")
+        ).expec()
+
+    def get_heuristic_variables(self) -> List[Variable]:
+        return [
+            float_variable(
+                "distance_between_target_and_generating",
+                lower_bound=literal(0),
+                structure=TIME_AND_SCENARIO_FREE,
+            ),
+            float_variable(
+                "max_distance_between_target_and_generating",
+                lower_bound=literal(0),
+                structure=CONSTANT_PER_SCENARIO,
+            ),
+            float_variable(
+                "violation_lower_rule_curve",
+                lower_bound=literal(0),
+                structure=TIME_AND_SCENARIO_FREE,
+            ),
+            float_variable(
+                "violation_upper_rule_curve",
+                lower_bound=literal(0),
+                structure=TIME_AND_SCENARIO_FREE,
+            ),
+            float_variable(
+                "max_violation_lower_rule_curve",
+                lower_bound=literal(0),
+                structure=CONSTANT_PER_SCENARIO,
+            ),
+            float_variable(
+                "gap_to_target",
+                lower_bound=literal(0),
+                structure=CONSTANT_PER_SCENARIO,
+            ),
+        ]
+
+    def get_heuristic_parameters(self) -> List[Parameter]:
+        return [
+            float_parameter("generating_target", TIME_AND_SCENARIO_FREE),
+            float_parameter("overall_target", CONSTANT_PER_SCENARIO),
+        ]
+
+    def get_heuristic_constraints(self) -> List[Constraint]:
+        list_constraint = [
+            Constraint(
+                "Respect generating target",
+                var("generating").sum() + var("gap_to_target")
+                == param("overall_target"),
+            ),
+            Constraint(
+                "Definition of distance between target and generating",
+                var("distance_between_target_and_generating")
+                >= param("generating_target") - var("generating"),
+            ),
+            Constraint(
+                "Definition of distance between generating and target",
+                var("distance_between_target_and_generating")
+                >= var("generating") - param("generating_target"),
+            ),
+            Constraint(
+                "Definition of max distance between generating and target",
+                var("max_distance_between_target_and_generating")
+                >= var("distance_between_target_and_generating"),
+            ),
+            Constraint(
+                "Definition of violation of lower rule curve",
+                var("violation_lower_rule_curve") + var("level")
+                >= param("lower_rule_curve"),
+            ),
+            Constraint(
+                "Definition of violation of upper rule curve",
+                var("violation_upper_rule_curve") - var("level")
+                >= -param("upper_rule_curve"),
+            ),
+            Constraint(
+                "Definition of max violation of lower rule curve",
+                var("max_violation_lower_rule_curve")
+                >= var("violation_lower_rule_curve"),
+            ),
+        ]
+        if self.horizon == "monthly":
+            list_constraint.append(
+                Constraint("No overflow", var("overflow") <= literal(0))
+            )
+            list_constraint.append(
+                Constraint("No gap to target", var("gap_to_target") <= literal(0))
+            )
+
+        return list_constraint