diff --git a/src/calculate_reward_and_bellman_values.py b/src/calculate_reward_and_bellman_values.py index a62ded3..9720868 100644 --- a/src/calculate_reward_and_bellman_values.py +++ b/src/calculate_reward_and_bellman_values.py @@ -1,219 +1,11 @@ from itertools import product import numpy as np -from scipy.interpolate import interp1d from estimation import PieceWiseLinearInterpolator, RewardApproximation -from read_antares_data import Reservoir, TimeScenarioIndex, TimeScenarioParameter -from type_definition import Array1D, Array2D, Callable, Dict, Iterable, List, Optional - - -class ReservoirManagement: - - def __init__( - self, - reservoir: Reservoir, - penalty_bottom_rule_curve: float = 0, - penalty_upper_rule_curve: float = 0, - penalty_final_level: float = 0, - force_final_level: bool = False, - final_level: Optional[float] = None, - overflow: bool = True, - ) -> None: - """Class to describe reservoir management parameters. - - Args: - reservoir (Reservoir): Reservoir in question - penalty_bottom_rule_curve (float, optional): Penalty for violating bottom rule curve. Defaults to 0. - penalty_upper_rule_curve (float, optional): Penalty for violating upper rule curve. Defaults to 0. - penalty_final_level (float, optional): Penalty for not respecting final level. Defaults to 0. - force_final_level (bool, optional): Whether final level is imposed. Defaults to False. - final_level (Optional[float], optional): Final level to impose, if not specified is equal to initial level. Defaults to None. - overflow (bool, optional) : Whether overflow is possible or forbiden. Defaults to True. - """ - - self.reservoir = reservoir - self.penalty_bottom_rule_curve = penalty_bottom_rule_curve - self.penalty_upper_rule_curve = penalty_upper_rule_curve - self.overflow = overflow - - if force_final_level: - self.penalty_final_level = penalty_final_level - if final_level: - self.final_level = final_level - else: - self.final_level = reservoir.initial_level - else: - self.final_level = False - - def get_penalty(self, week: int, len_week: int) -> Callable: - """ - Return a function to evaluate penalities for violating rule curves for any level of stock. - - Parameters - ---------- - week:int : - Week considered - len_week:int : - Total number of weeks - - Returns - ------- - - """ - if week == len_week - 1 and self.final_level: - pen = interp1d( - [ - 0, - self.final_level, - self.reservoir.capacity, - ], - [ - -self.penalty_final_level * (self.final_level), - 0, - -self.penalty_final_level - * (self.reservoir.capacity - self.final_level), - ], - ) - else: - pen = interp1d( - [ - 0, - self.reservoir.bottom_rule_curve[week], - self.reservoir.upper_rule_curve[week], - self.reservoir.capacity, - ], - [ - -self.penalty_bottom_rule_curve - * (self.reservoir.bottom_rule_curve[week]), - 0, - 0, - -self.penalty_upper_rule_curve - * (self.reservoir.capacity - self.reservoir.upper_rule_curve[week]), - ], - ) - return pen - - -class MultiStockManagement: - - def __init__(self, list_reservoirs: List[ReservoirManagement]) -> None: - """Describes reservoir management for all stocks - - Args: - list_reservoirs (List[ReservoirManagement]): List of reservoir management - """ - self.dict_reservoirs = {} - for res in list_reservoirs: - self.dict_reservoirs[res.reservoir.area] = res - - def get_disc( - self, - method: str, - week: int, - xNsteps: int, - reference_pt: np.ndarray, - correlation_matrix: np.ndarray, - alpha: float = 1.0, - in_out_ratio: float = 3.0, - ) -> Array1D: - n_reservoirs = len(self.dict_reservoirs) - if len(reference_pt.shape) > 1: - reference_pt = np.mean(reference_pt, axis=1) - lbs = np.array( - [ - mng.reservoir.bottom_rule_curve[week] * alpha - for mng in self.dict_reservoirs.values() - ] - ) - ubs = np.array( - [ - mng.reservoir.upper_rule_curve[week] * alpha - + mng.reservoir.capacity * (1 - alpha) - for mng in self.dict_reservoirs.values() - ] - ) - full = np.array( - [mng.reservoir.capacity for mng in self.dict_reservoirs.values()] - ) - empty = np.array([0] * n_reservoirs) - n_pts_above = np.maximum( - 1 + 2 * (full - ubs > 0), - np.round( - xNsteps - * (full - ubs) - / (full - empty + (in_out_ratio - 1) * (ubs - lbs)) - ).astype(int), - ) - n_pts_in = np.round( - xNsteps - * (ubs - lbs) - * in_out_ratio - / (full - empty + (in_out_ratio - 1) * (ubs - lbs)) - ).astype(int) - n_pts_below = np.maximum( - 1, - np.round( - xNsteps - * (lbs - empty) - / (full - empty + (in_out_ratio - 1) * (ubs - lbs)) - ).astype(int), - ) - n_pts_in += xNsteps - ( - n_pts_below + n_pts_in + n_pts_above - ) # Make sure total adds up - if method == "lines": - above_curve_pts = [ - np.linspace(ubs[r], full[r], n_pts_above[r], endpoint=True) - for r in range(n_reservoirs) - ] - in_curve_pts = [ - np.linspace(lbs[r], ubs[r], n_pts_in[r], endpoint=False) - for r in range(n_reservoirs) - ] - below_curve_pts = [ - np.linspace(empty[r], lbs[r], n_pts_below[r], endpoint=False) - for r in range(n_reservoirs) - ] - all_pts = np.array( - [ - np.concatenate( - (below_curve_pts[r], in_curve_pts[r], above_curve_pts[r]) - ) - for r in range(n_reservoirs) - ] - ).T - diffs_to_ref = all_pts[:, None] - reference_pt[None, :] # Disc * R - diffs_to_ref = ( - diffs_to_ref[:, :, None] * np.eye(n_reservoirs)[None, :, :] - ) # Disc * R * R - diffs_to_ref = np.dot(diffs_to_ref, correlation_matrix) # Disc * R * R - new_levels = reference_pt[None, None, :] + diffs_to_ref # Disc * R * R - new_levels = np.maximum( - new_levels, empty[None, None, :] - ) # Do or do not make other points leave guiding curves ? - new_levels = np.minimum( - new_levels, full[None, None, :] - ) # We chose the former - levels = np.reshape( - new_levels, (xNsteps * n_reservoirs, n_reservoirs) - ) # (Disc * R) * R - else: - # Listing all levels - levels_discretization = product( - *[ - np.concatenate( - [ - [0], - np.linspace(lbs[i], ubs[i], xNsteps - 2), - [(alpha + 1) / 2 * mng.reservoir.capacity], - ] - ) - for i, mng in enumerate(self.dict_reservoirs.values()) - ] - ) - levels = np.array([level for level in levels_discretization]) - return levels +from read_antares_data import TimeScenarioIndex, TimeScenarioParameter +from reservoir_management import ReservoirManagement +from type_definition import Array1D, Callable, Dict, Iterable, List class BellmanValueCalculation: diff --git a/src/display.py b/src/display.py index 4de5600..22d879b 100644 --- a/src/display.py +++ b/src/display.py @@ -2,8 +2,8 @@ import plotly.graph_objects as go from tqdm import tqdm -from multi_stock_bellman_value_calculation import MultiStockManagement from read_antares_data import TimeScenarioParameter +from reservoir_management import MultiStockManagement from type_definition import Array1D, Dict diff --git a/src/estimation.py b/src/estimation.py index 96d7542..fab7802 100644 --- a/src/estimation.py +++ b/src/estimation.py @@ -4,6 +4,8 @@ from hyperplane_decomposition import decompose_hyperplanes from hyperplane_interpolation import get_interpolation from read_antares_data import TimeScenarioIndex, TimeScenarioParameter +from reservoir_management import MultiStockManagement +from stock_discretization import StockDiscretization from type_definition import Array1D, Callable, Dict, List, Optional, Union @@ -12,6 +14,19 @@ class Estimator: def __init__(self) -> None: raise NotImplemented + def update( + self, + Vu: float, + slope: Dict[str, float], + n_scenario: int, + idx: int, + list_areas: List[str], + ) -> None: + raise NotImplemented + + def get_value(self, x: Dict[str, float]) -> float: + return NotImplemented + class PieceWiseLinearInterpolator: @@ -39,15 +54,64 @@ def __init__( def __getitem__(self, key: str) -> PieceWiseLinearInterpolator: return self.estimators[key] + def update( + self, + Vu: float, + slope: Dict[str, float], + n_scenario: int, + idx: int, + list_areas: List[str], + ) -> None: + for area in list_areas: + self.estimators[area].costs[idx] += -Vu / n_scenario + + def get_value(self, x: Dict[str, float]) -> float: + return -sum([self.estimators[area](y) for area, y in x.items()]) + class BellmanValueEstimation(Estimator): - def __init__(self, value: Dict[str, Array1D]): + def __init__( + self, value: Dict[str, Array1D], stock_discretization: StockDiscretization + ): self.V = value + self.discretization = stock_discretization def __getitem__(self, key: str) -> Array1D: return self.V[key] + def update( + self, + Vu: float, + slope: Dict[str, float], + n_scenario: int, + idx: int, + list_areas: List[str], + ) -> None: + self.V["intercept"][idx] += Vu / n_scenario + for area in list_areas: + self.V[f"slope_{area}"][idx] += slope[area] / n_scenario + + def get_value(self, x: Dict[str, float]) -> float: + return max( + [ + sum( + [ + self.V[f"slope_{area}"][idx] + * ( + x[area] + - self.discretization.list_discretization[area][idx[i]] + ) + for i, area in enumerate( + self.discretization.list_discretization.keys() + ) + ] + ) + + self.V["intercept"][idx] + for idx in self.discretization.get_product_stock_discretization() + ] + ) + class MultiVariateEstimator: diff --git a/src/functions_iterative.py b/src/functions_iterative.py index 49e4cb3..a4e33a3 100644 --- a/src/functions_iterative.py +++ b/src/functions_iterative.py @@ -6,17 +6,18 @@ from calculate_reward_and_bellman_values import ( BellmanValueCalculation, MultiStockBellmanValueCalculation, - MultiStockManagement, ReservoirManagement, RewardApproximation, ) from estimation import ( BellmanValueEstimation, + Estimator, PieceWiseLinearInterpolator, UniVariateEstimator, ) from optimization import AntaresProblem, Basis from read_antares_data import TimeScenarioIndex, TimeScenarioParameter +from reservoir_management import MultiStockManagement from type_definition import Array1D, Array2D, Array3D, Array4D, Dict, List, Union @@ -89,7 +90,7 @@ def compute_upper_bound( param: TimeScenarioParameter, multi_bellman_value_calculation: MultiStockBellmanValueCalculation, list_models: Dict[TimeScenarioIndex, AntaresProblem], - V: Dict[int, Union[UniVariateEstimator, BellmanValueEstimation]], + V: Dict[int, Estimator], ) -> tuple[float, Dict[TimeScenarioIndex, Dict[str, float]], Array3D]: """ Compute an approximate upper bound on the overall problem by solving the real complete Antares problem with Bellman values. diff --git a/src/launch_calculation.py b/src/launch_calculation.py index afd49c6..901831b 100644 --- a/src/launch_calculation.py +++ b/src/launch_calculation.py @@ -1,4 +1,7 @@ +import numpy as np + from calculate_reward_and_bellman_values import ReservoirManagement +from estimation import UniVariateEstimator from functions_iterative import itr_control from multi_stock_bellman_value_calculation import MultiStockManagement from read_antares_data import TimeScenarioParameter @@ -41,13 +44,26 @@ def calculate_bellman_values( if method == "direct": # Compute Bellman values directly - vb = calculate_bellman_value_directly( + intermediate_vb, _, _ = calculate_bellman_value_directly( param=param, multi_stock_management=MultiStockManagement([reservoir_management]), output_path=output_path, X={reservoir_management.reservoir.area: X}, solver=solver, - )[reservoir_management.reservoir.area] + univariate=True, + ) + + vb = np.transpose( + [ + [ + intermediate_vb[week].get_value( + {reservoir_management.reservoir.area: x} + ) + for x in X + ] + for week in range(param.len_week + 1) + ] + ) elif method == "precalculated": # or with precalulated reward diff --git a/src/multi_stock_bellman_value_calculation.py b/src/multi_stock_bellman_value_calculation.py index ebdc8fa..cd13fe6 100644 --- a/src/multi_stock_bellman_value_calculation.py +++ b/src/multi_stock_bellman_value_calculation.py @@ -15,7 +15,6 @@ from calculate_reward_and_bellman_values import ( BellmanValueCalculation, MultiStockBellmanValueCalculation, - MultiStockManagement, RewardApproximation, ) from display import ConvergenceProgressBar, draw_usage_values, draw_uvs_sddp @@ -28,6 +27,8 @@ solve_for_optimal_trajectory, ) from read_antares_data import TimeScenarioIndex, TimeScenarioParameter +from reservoir_management import MultiStockManagement +from stock_discretization import StockDiscretization jl = juliacall.Main @@ -35,139 +36,6 @@ Array2D = Annotated[npt.NDArray[np.float32], Literal["N", "N"]] -def calculate_bellman_value_multi_stock( - param: TimeScenarioParameter, - multi_stock_management: MultiStockManagement, - output_path: str, - X: Dict[str, Array1D], - name_solver: str = "CLP", -) -> tuple[Dict[int, Dict[str, npt.NDArray[np.float32]]], float, float]: - """Function to calculate multivariate Bellman values - - Args: - param (TimeScenarioParameter): Time-related parameters for the Antares study - multi_stock_management (MultiStockManagement): Stocks considered for Bellman values - output_path (str): Path to mps files describing optimization problems - X (Dict[str, Array1D]): Discretization of sotck levels for Bellman values for each reservoir - name_solver (str, optional): Solver to use with ortools. Defaults to "CLP". - - Returns: - Dict[int, Dict[str, npt.NDArray[np.float32]]]: Bellman values for each week. Bellman values are described by a slope for each area and a intercept - """ - - list_models: Dict[TimeScenarioIndex, AntaresProblem] = {} - for week in range(param.len_week): - for scenario in range(param.len_scenario): - m = AntaresProblem( - scenario=scenario, - week=week, - path=output_path, - itr=1, - name_solver=name_solver, - name_scenario=param.name_scenario[scenario], - ) - - m.create_weekly_problem_itr( - param=param, - multi_stock_management=multi_stock_management, - ) - list_models[TimeScenarioIndex(week, scenario)] = m - - reward: Dict[str, Dict[TimeScenarioIndex, RewardApproximation]] = {} - for area, reservoir_management in multi_stock_management.dict_reservoirs.items(): - reward[area] = {} - for week in range(param.len_week): - for scenario in range(param.len_scenario): - r = RewardApproximation( - lb_control=-reservoir_management.reservoir.max_pumping[week], - ub_control=reservoir_management.reservoir.max_generating[week], - ub_reward=0, - ) - reward[area][TimeScenarioIndex(week, scenario)] = r - - bellman_value_calculation = [] - for area, reservoir_management in multi_stock_management.dict_reservoirs.items(): - bellman_value_calculation.append( - BellmanValueCalculation( - param=param, - reward=reward[area], - reservoir_management=reservoir_management, - stock_discretization=X[area], - ) - ) - multi_bellman_value_calculation = MultiStockBellmanValueCalculation( - bellman_value_calculation - ) - - dim_bellman_value = tuple(len(x) for x in X.values()) - V = { - week: { - name: np.zeros((dim_bellman_value), dtype=np.float32) - for name in ["intercept"] - + [ - f"slope_{area}" - for area in multi_bellman_value_calculation.dict_reservoirs.keys() - ] - } - for week in range(param.len_week + 1) - } - - for week in range(param.len_week - 1, -1, -1): - for scenario in range(param.len_scenario): - print(f"{week} {scenario}", end="\r") - m = list_models[TimeScenarioIndex(week, scenario)] - - for ( - idx - ) in multi_bellman_value_calculation.get_product_stock_discretization(): - _, _, Vu, slope, _, _, _ = m.solve_problem_with_bellman_values( - multi_bellman_value_calculation=multi_bellman_value_calculation, - V=BellmanValueEstimation(V[week + 1]), - level_i={ - area: multi_bellman_value_calculation.dict_reservoirs[ - area - ].stock_discretization[idx[i]] - for i, area in enumerate(m.range_reservoir) - }, - take_into_account_z_and_y=True, - ) - V[week]["intercept"][idx] += Vu / param.len_scenario - for area in m.range_reservoir: - V[week][f"slope_{area}"][idx] += slope[area] / param.len_scenario - - lower_bound = max( - [ - sum( - [ - V[0][f"slope_{area}"][idx] - * ( - multi_bellman_value_calculation.dict_reservoirs[ - area - ].reservoir_management.reservoir.initial_level - - multi_bellman_value_calculation.dict_reservoirs[ - area - ].stock_discretization[idx[i]] - ) - for i, area in enumerate( - multi_bellman_value_calculation.dict_reservoirs.keys() - ) - ] - ) - + V[0]["intercept"][idx] - for idx in multi_bellman_value_calculation.get_product_stock_discretization() - ] - ) - - upper_bound, _, _ = compute_upper_bound( - param=param, - multi_bellman_value_calculation=multi_bellman_value_calculation, - list_models=list_models, - V={week: BellmanValueEstimation(V[week]) for week in range(param.len_week + 1)}, - ) - - return V, lower_bound, upper_bound - - def initialize_antares_problems( param: TimeScenarioParameter, multi_stock_management: MultiStockManagement, diff --git a/src/optimization.py b/src/optimization.py index 6d88d49..ff26b44 100644 --- a/src/optimization.py +++ b/src/optimization.py @@ -6,22 +6,18 @@ import numpy as np import ortools.linear_solver.pywraplp as pywraplp from ortools.linear_solver.python import model_builder -from scipy.interpolate import interp1d -from calculate_reward_and_bellman_values import ( - BellmanValueCalculation, - MultiStockBellmanValueCalculation, - MultiStockManagement, -) +from calculate_reward_and_bellman_values import MultiStockBellmanValueCalculation from estimation import ( BellmanValueEstimation, + Estimator, LinearCostEstimator, LinearInterpolator, - PieceWiseLinearInterpolator, UniVariateEstimator, ) from read_antares_data import TimeScenarioIndex, TimeScenarioParameter -from type_definition import Array1D, Array2D, Dict, List, Union, npt +from reservoir_management import MultiStockManagement +from type_definition import Array1D, Dict, List, Union class Basis: @@ -469,7 +465,7 @@ def set_constraints_predefined_control(self, control: float, area: str) -> None: def set_constraints_initial_level_and_bellman_values( self, - V: Union[UniVariateEstimator, BellmanValueEstimation], + V: Estimator, all_level_i: Dict[str, float], multi_bellman_value_calculation: MultiStockBellmanValueCalculation, ) -> List[pywraplp.Constraint]: @@ -652,7 +648,7 @@ def get_solution_value(self, name_variable: str) -> Dict[str, float]: def solve_problem_with_bellman_values( self, multi_bellman_value_calculation: MultiStockBellmanValueCalculation, - V: Union[UniVariateEstimator, BellmanValueEstimation], + V: Estimator, level_i: Dict[str, float], take_into_account_z_and_y: bool, find_optimal_basis: bool = True, diff --git a/src/reservoir_management.py b/src/reservoir_management.py new file mode 100644 index 0000000..49da810 --- /dev/null +++ b/src/reservoir_management.py @@ -0,0 +1,215 @@ +from itertools import product + +import numpy as np +from scipy.interpolate import interp1d + +from read_antares_data import Reservoir +from type_definition import Array1D, Callable, List, Optional + + +class ReservoirManagement: + + def __init__( + self, + reservoir: Reservoir, + penalty_bottom_rule_curve: float = 0, + penalty_upper_rule_curve: float = 0, + penalty_final_level: float = 0, + force_final_level: bool = False, + final_level: Optional[float] = None, + overflow: bool = True, + ) -> None: + """Class to describe reservoir management parameters. + + Args: + reservoir (Reservoir): Reservoir in question + penalty_bottom_rule_curve (float, optional): Penalty for violating bottom rule curve. Defaults to 0. + penalty_upper_rule_curve (float, optional): Penalty for violating upper rule curve. Defaults to 0. + penalty_final_level (float, optional): Penalty for not respecting final level. Defaults to 0. + force_final_level (bool, optional): Whether final level is imposed. Defaults to False. + final_level (Optional[float], optional): Final level to impose, if not specified is equal to initial level. Defaults to None. + overflow (bool, optional) : Whether overflow is possible or forbiden. Defaults to True. + """ + + self.reservoir = reservoir + self.penalty_bottom_rule_curve = penalty_bottom_rule_curve + self.penalty_upper_rule_curve = penalty_upper_rule_curve + self.overflow = overflow + + if force_final_level: + self.penalty_final_level = penalty_final_level + if final_level: + self.final_level = final_level + else: + self.final_level = reservoir.initial_level + else: + self.final_level = False + + def get_penalty(self, week: int, len_week: int) -> Callable: + """ + Return a function to evaluate penalities for violating rule curves for any level of stock. + + Parameters + ---------- + week:int : + Week considered + len_week:int : + Total number of weeks + + Returns + ------- + + """ + if week == len_week - 1 and self.final_level: + pen = interp1d( + [ + 0, + self.final_level, + self.reservoir.capacity, + ], + [ + -self.penalty_final_level * (self.final_level), + 0, + -self.penalty_final_level + * (self.reservoir.capacity - self.final_level), + ], + ) + else: + pen = interp1d( + [ + 0, + self.reservoir.bottom_rule_curve[week], + self.reservoir.upper_rule_curve[week], + self.reservoir.capacity, + ], + [ + -self.penalty_bottom_rule_curve + * (self.reservoir.bottom_rule_curve[week]), + 0, + 0, + -self.penalty_upper_rule_curve + * (self.reservoir.capacity - self.reservoir.upper_rule_curve[week]), + ], + ) + return pen + + +class MultiStockManagement: + + def __init__(self, list_reservoirs: List[ReservoirManagement]) -> None: + """Describes reservoir management for all stocks + + Args: + list_reservoirs (List[ReservoirManagement]): List of reservoir management + """ + self.dict_reservoirs = {} + for res in list_reservoirs: + self.dict_reservoirs[res.reservoir.area] = res + + def get_disc( + self, + method: str, + week: int, + xNsteps: int, + reference_pt: np.ndarray, + correlation_matrix: np.ndarray, + alpha: float = 1.0, + in_out_ratio: float = 3.0, + ) -> Array1D: + n_reservoirs = len(self.dict_reservoirs) + if len(reference_pt.shape) > 1: + reference_pt = np.mean(reference_pt, axis=1) + lbs = np.array( + [ + mng.reservoir.bottom_rule_curve[week] * alpha + for mng in self.dict_reservoirs.values() + ] + ) + ubs = np.array( + [ + mng.reservoir.upper_rule_curve[week] * alpha + + mng.reservoir.capacity * (1 - alpha) + for mng in self.dict_reservoirs.values() + ] + ) + full = np.array( + [mng.reservoir.capacity for mng in self.dict_reservoirs.values()] + ) + empty = np.array([0] * n_reservoirs) + n_pts_above = np.maximum( + 1 + 2 * (full - ubs > 0), + np.round( + xNsteps + * (full - ubs) + / (full - empty + (in_out_ratio - 1) * (ubs - lbs)) + ).astype(int), + ) + n_pts_in = np.round( + xNsteps + * (ubs - lbs) + * in_out_ratio + / (full - empty + (in_out_ratio - 1) * (ubs - lbs)) + ).astype(int) + n_pts_below = np.maximum( + 1, + np.round( + xNsteps + * (lbs - empty) + / (full - empty + (in_out_ratio - 1) * (ubs - lbs)) + ).astype(int), + ) + n_pts_in += xNsteps - ( + n_pts_below + n_pts_in + n_pts_above + ) # Make sure total adds up + if method == "lines": + above_curve_pts = [ + np.linspace(ubs[r], full[r], n_pts_above[r], endpoint=True) + for r in range(n_reservoirs) + ] + in_curve_pts = [ + np.linspace(lbs[r], ubs[r], n_pts_in[r], endpoint=False) + for r in range(n_reservoirs) + ] + below_curve_pts = [ + np.linspace(empty[r], lbs[r], n_pts_below[r], endpoint=False) + for r in range(n_reservoirs) + ] + all_pts = np.array( + [ + np.concatenate( + (below_curve_pts[r], in_curve_pts[r], above_curve_pts[r]) + ) + for r in range(n_reservoirs) + ] + ).T + diffs_to_ref = all_pts[:, None] - reference_pt[None, :] # Disc * R + diffs_to_ref = ( + diffs_to_ref[:, :, None] * np.eye(n_reservoirs)[None, :, :] + ) # Disc * R * R + diffs_to_ref = np.dot(diffs_to_ref, correlation_matrix) # Disc * R * R + new_levels = reference_pt[None, None, :] + diffs_to_ref # Disc * R * R + new_levels = np.maximum( + new_levels, empty[None, None, :] + ) # Do or do not make other points leave guiding curves ? + new_levels = np.minimum( + new_levels, full[None, None, :] + ) # We chose the former + levels = np.reshape( + new_levels, (xNsteps * n_reservoirs, n_reservoirs) + ) # (Disc * R) * R + else: + # Listing all levels + levels_discretization = product( + *[ + np.concatenate( + [ + [0], + np.linspace(lbs[i], ubs[i], xNsteps - 2), + [(alpha + 1) / 2 * mng.reservoir.capacity], + ] + ) + for i, mng in enumerate(self.dict_reservoirs.values()) + ] + ) + levels = np.array([level for level in levels_discretization]) + return levels diff --git a/src/simple_bellman_value_calculation.py b/src/simple_bellman_value_calculation.py index 9991e53..d9e1923 100644 --- a/src/simple_bellman_value_calculation.py +++ b/src/simple_bellman_value_calculation.py @@ -1,18 +1,23 @@ import numpy as np -from scipy.interpolate import interp1d from calculate_reward_and_bellman_values import ( BellmanValueCalculation, MultiStockBellmanValueCalculation, - MultiStockManagement, ReservoirManagement, RewardApproximation, ) -from estimation import PieceWiseLinearInterpolator, UniVariateEstimator +from estimation import ( + BellmanValueEstimation, + Estimator, + PieceWiseLinearInterpolator, + UniVariateEstimator, +) from functions_iterative import compute_upper_bound from optimization import AntaresProblem, Basis from read_antares_data import TimeScenarioIndex, TimeScenarioParameter -from type_definition import Array1D, Array2D, Dict +from reservoir_management import MultiStockManagement +from stock_discretization import StockDiscretization +from type_definition import Array1D, Array2D, Dict, Union def calculate_complete_reward( @@ -166,8 +171,9 @@ def calculate_bellman_value_directly( multi_stock_management: MultiStockManagement, output_path: str, X: Dict[str, Array1D], + univariate: bool, solver: str = "CLP", -) -> Dict[str, Array2D]: +) -> tuple[Dict[int, Estimator], float, float]: """ Algorithm to evaluate Bellman values directly. @@ -235,16 +241,37 @@ def calculate_bellman_value_directly( multi_bellman_value_calculation = MultiStockBellmanValueCalculation( bellman_value_calculation ) - - V = { - area: { - week: PieceWiseLinearInterpolator( - X[area], np.zeros(len(X[area]), dtype=np.float32) + V: Dict[int, Estimator] = {} + if univariate: + assert len(X) == 1 + dim_bellman_value = tuple(len(x) for x in X.values()) + V = { + week: UniVariateEstimator( + { + area: PieceWiseLinearInterpolator( + X[area], np.zeros(len(X[area]), dtype=np.float32) + ) + for area in multi_stock_management.dict_reservoirs.keys() + } + ) + for week in range(param.len_week + 1) + } + else: + dim_bellman_value = tuple(len(x) for x in X.values()) + V = { + week: BellmanValueEstimation( + { + name: np.zeros((dim_bellman_value), dtype=np.float32) + for name in ["intercept"] + + [ + f"slope_{area}" + for area in multi_bellman_value_calculation.dict_reservoirs.keys() + ] + }, + StockDiscretization(X), ) for week in range(param.len_week + 1) } - for area in multi_stock_management.dict_reservoirs.keys() - } for week in range(param.len_week - 1, -1, -1): for scenario in range(param.len_scenario): @@ -254,11 +281,9 @@ def calculate_bellman_value_directly( for ( idx ) in multi_bellman_value_calculation.get_product_stock_discretization(): - _, _, Vu, _, _, _, _ = m.solve_problem_with_bellman_values( + _, _, Vu, slope, _, _, _ = m.solve_problem_with_bellman_values( multi_bellman_value_calculation=multi_bellman_value_calculation, - V=UniVariateEstimator( - {area: V[area][week + 1] for area in m.range_reservoir} - ), + V=V[week + 1], level_i={ area: multi_bellman_value_calculation.dict_reservoirs[ area @@ -268,36 +293,29 @@ def calculate_bellman_value_directly( find_optimal_basis=False, take_into_account_z_and_y=True, ) - V[area][week].costs[idx] += -Vu / param.len_scenario + V[week].update( + Vu, + slope, + param.len_scenario, + idx, + list(multi_stock_management.dict_reservoirs.keys()), + ) - V0 = sum( - [ - V[area][0]( - multi_stock_management.dict_reservoirs[area].reservoir.initial_level - ) + V0 = V[0].get_value( + { + area: multi_stock_management.dict_reservoirs[area].reservoir.initial_level for area in multi_stock_management.dict_reservoirs.keys() - ] + } ) upper_bound, controls, current_itr = compute_upper_bound( param=param, multi_bellman_value_calculation=multi_bellman_value_calculation, list_models=list_models, - V={ - week: UniVariateEstimator( - { - area: V[area][week] - for area in multi_stock_management.dict_reservoirs.keys() - } - ) - for week in range(param.len_week + 1) - }, + V={week: V[week] for week in range(param.len_week + 1)}, ) gap = upper_bound + V0 print(gap, upper_bound, -V0) - return { - area: np.transpose([V[area][week].costs for week in range(param.len_week + 1)]) - for area in multi_stock_management.dict_reservoirs.keys() - } + return {week: V[week] for week in range(param.len_week + 1)}, V0, upper_bound diff --git a/src/stock_discretization.py b/src/stock_discretization.py new file mode 100644 index 0000000..a406f41 --- /dev/null +++ b/src/stock_discretization.py @@ -0,0 +1,14 @@ +from itertools import product + +from type_definition import Array1D, Dict, Iterable + + +class StockDiscretization: + + def __init__(self, list_discretization: Dict[str, Array1D]) -> None: + self.list_discretization = list_discretization + + def get_product_stock_discretization(self) -> Iterable: + return product( + *[[i for i in range(len(x))] for x in self.list_discretization.values()] + ) diff --git a/tests/test_basis.py b/tests/test_basis.py index e534686..37f19e6 100644 --- a/tests/test_basis.py +++ b/tests/test_basis.py @@ -2,7 +2,6 @@ import ortools.linear_solver.pywraplp as pywraplp import pytest -from calculate_reward_and_bellman_values import MultiStockManagement from estimation import PieceWiseLinearInterpolator, UniVariateEstimator from functions_iterative import ( BellmanValueCalculation, @@ -15,6 +14,7 @@ from multi_stock_bellman_value_calculation import MultiStockBellmanValueCalculation from optimization import AntaresProblem, Basis from read_antares_data import Reservoir +from reservoir_management import MultiStockManagement def test_basis_with_xpress() -> None: diff --git a/tests/test_bellman_value_exact.py b/tests/test_bellman_value_exact.py index fb03fa8..2c97c3d 100644 --- a/tests/test_bellman_value_exact.py +++ b/tests/test_bellman_value_exact.py @@ -1,15 +1,13 @@ import numpy as np import pytest +from estimation import BellmanValueEstimation from functions_iterative import ReservoirManagement, TimeScenarioParameter -from multi_stock_bellman_value_calculation import ( - MultiStockManagement, - calculate_bellman_value_multi_stock, -) +from multi_stock_bellman_value_calculation import MultiStockManagement from read_antares_data import Reservoir from simple_bellman_value_calculation import calculate_bellman_value_directly -expected_vb = np.array( +expected_vb = -np.array( [ [ -5.8881910e09, @@ -174,6 +172,171 @@ ], ) +expected_vb_lower_approximation = np.array( + [ + [ + 5.88819046e09, + 5.37158298e09, + 4.30354534e09, + 3.62174925e09, + 1.99857485e09, + 0.00000000e00, + ], + [ + 5.28486810e09, + 4.36279347e09, + 3.38627744e09, + 2.42429722e09, + 1.40291200e09, + 0.00000000e00, + ], + [ + 5.12697223e09, + 4.20489754e09, + 3.22838170e09, + 2.26640128e09, + 1.24501632e09, + 0.00000000e00, + ], + [ + 4.98446182e09, + 4.04738202e09, + 3.07425285e09, + 2.10850586e09, + 1.10278950e09, + 0.00000000e00, + ], + [ + 4.87915648e09, + 3.94207644e09, + 2.96894746e09, + 1.99659827e09, + 9.97484032e08, + 0.00000000e00, + ], + [ + 4.77385114e09, + 3.83677107e09, + 2.86364211e09, + 1.89129290e09, + 8.92178688e08, + 0.00000000e00, + ], + [ + 4.66854579e09, + 3.73146573e09, + 2.75833680e09, + 1.78598758e09, + 7.86873472e08, + 0.00000000e00, + ], + [ + 4.56324045e09, + 3.62616042e09, + 2.65303149e09, + 1.68068227e09, + 6.97788608e08, + 0.00000000e00, + ], + [ + 4.45793514e09, + 3.52085510e09, + 2.54772621e09, + 1.57537709e09, + 6.45156815e08, + 0.00000000e00, + ], + [ + 4.35262982e09, + 3.41554982e09, + 2.44242106e09, + 1.47007194e09, + 5.92525120e08, + 0.00000000e00, + ], + [ + 4.24732451e09, + 3.31024467e09, + 2.33711593e09, + 1.39639194e09, + 5.39893440e08, + 0.00000000e00, + ], + [ + 4.14201938e09, + 3.20493954e09, + 2.23181082e09, + 1.34376014e09, + 4.87261792e08, + 0.00000000e00, + ], + [ + 4.03671424e09, + 3.09963443e09, + 2.13425011e09, + 1.29112834e09, + 4.34630153e08, + 0.00000000e00, + ], + [ + 3.93140910e09, + 2.99432932e09, + 2.08161830e09, + 1.23849656e09, + 3.81998560e08, + 0.00000000e00, + ], + [ + 3.82610397e09, + 2.88902421e09, + 2.02898650e09, + 1.18586485e09, + 3.29367104e08, + 0.00000000e00, + ], + [ + 3.75680376e09, + 2.81986609e09, + 1.97635470e09, + 1.13323315e09, + 2.76736128e08, + 0.00000000e00, + ], + [ + 3.70417203e09, + 2.76723430e09, + 1.92372291e09, + 1.08060146e09, + 2.24105184e08, + 0.00000000e00, + ], + [ + 3.65154032e09, + 2.71460252e09, + 1.87109120e09, + 1.02796976e09, + 1.71474288e08, + 0.00000000e00, + ], + [ + 3.59890867e09, + 2.66197074e09, + 1.81845952e09, + 9.75338066e08, + 1.18843432e08, + 0.00000000e00, + ], + [ + 3.54627702e09, + 2.60933898e09, + 1.76582784e09, + 9.22706432e08, + 6.62126120e07, + 0.00000000e00, + ], + ] +) + def test_bellman_value_exact() -> None: @@ -189,14 +352,24 @@ def test_bellman_value_exact() -> None: xNsteps = 20 X = np.linspace(0, reservoir.capacity, num=xNsteps) - vb = calculate_bellman_value_directly( + vb, lb, ub = calculate_bellman_value_directly( param=param, multi_stock_management=MultiStockManagement([reservoir_management]), output_path="test_data/one_node", X={reservoir.area: X}, + univariate=True, ) - assert vb[reservoir.area] == pytest.approx(expected_vb) + assert lb == pytest.approx(4410021312) + + assert ub == pytest.approx(4410021125.56477) + + assert np.transpose( + [ + [vb[week].get_value({reservoir.area: x}) for x in X] + for week in range(param.len_week + 1) + ] + ) == pytest.approx(expected_vb) def test_bellman_value_exact_with_multi_stock() -> None: @@ -213,17 +386,26 @@ def test_bellman_value_exact_with_multi_stock() -> None: xNsteps = 20 X = np.linspace(0, reservoir.capacity, num=xNsteps) - vb, _, _ = calculate_bellman_value_multi_stock( + vb, lb, ub = calculate_bellman_value_directly( param=param, multi_stock_management=MultiStockManagement([reservoir_management]), output_path="test_data/one_node", X={"area": X}, + univariate=False, ) + assert lb == pytest.approx(4410021218.59294) + + assert ub == pytest.approx(4410021094.449415) + computed_vb = np.transpose( - [-vb[i]["intercept"] for i in range(param.len_week)] + [np.zeros(xNsteps)] + [ + [vb[week].get_value({reservoir.area: x}) for x in X] + for week in range(param.len_week + 1) + ] ) - # for week in range(param.len_week): - # assert computed_vb[:, param.len_week - week - 1] == pytest.approx( - # expected_vb[:, param.len_week - week - 1] - # ) + + for week in range(param.len_week): + assert computed_vb[:, param.len_week - week - 1] == pytest.approx( + expected_vb_lower_approximation[:, param.len_week - week - 1] + ) diff --git a/tests/test_bellman_value_exact_two_stocks.py b/tests/test_bellman_value_exact_two_stocks.py index 105bbcf..f4d0d03 100644 --- a/tests/test_bellman_value_exact_two_stocks.py +++ b/tests/test_bellman_value_exact_two_stocks.py @@ -6,7 +6,6 @@ from calculate_reward_and_bellman_values import ( BellmanValueCalculation, MultiStockBellmanValueCalculation, - MultiStockManagement, RewardApproximation, ) from estimation import BellmanValueEstimation @@ -15,12 +14,11 @@ TimeScenarioIndex, TimeScenarioParameter, ) -from multi_stock_bellman_value_calculation import ( - AntaresProblem, - Dict, - calculate_bellman_value_multi_stock, -) +from multi_stock_bellman_value_calculation import AntaresProblem, Dict from read_antares_data import Reservoir +from reservoir_management import MultiStockManagement +from simple_bellman_value_calculation import calculate_bellman_value_directly +from stock_discretization import StockDiscretization def test_iterate_over_stock_discretization() -> None: @@ -229,7 +227,7 @@ def test_solve_with_bellman_multi_stock() -> None: _, _, Vu, slope, _, xf, _ = m.solve_problem_with_bellman_values( multi_bellman_value_calculation=multi_bellman_value_calculation, - V=BellmanValueEstimation(V), + V=BellmanValueEstimation(V, StockDiscretization(X)), level_i={ area: multi_bellman_value_calculation.dict_reservoirs[ area @@ -279,16 +277,18 @@ def test_bellman_value_exact_calculation_multi_stock() -> None: x_1 = np.linspace(0, reservoir_1.capacity, num=5) x_2 = np.linspace(0, reservoir_2.capacity, num=5) - vb, lb, ub = calculate_bellman_value_multi_stock( + vb, lb, ub = calculate_bellman_value_directly( param=param, multi_stock_management=all_reservoirs, output_path="test_data/two_nodes", X={"area_1": x_1, "area_2": x_2}, + univariate=False, ) assert lb == pytest.approx(419088906.63159156) assert ub == pytest.approx(798837417.3288709) + assert type(vb[0]) is BellmanValueEstimation assert vb[0]["intercept"] == pytest.approx( np.array( [ @@ -376,3 +376,46 @@ def test_bellman_value_exact_calculation_multi_stock() -> None: ] ) ) + + +def test_bellman_value_exact_calculation_univariate() -> None: + + param = TimeScenarioParameter(len_week=5, len_scenario=1) + + reservoir_1 = Reservoir("test_data/two_nodes", "area_1") + reservoir_management_1 = ReservoirManagement( + reservoir=reservoir_1, + penalty_bottom_rule_curve=3000, + penalty_upper_rule_curve=3000, + penalty_final_level=3000, + force_final_level=True, + ) + + reservoir_2 = Reservoir("test_data/two_nodes", "area_2") + reservoir_management_2 = ReservoirManagement( + reservoir=reservoir_2, + penalty_bottom_rule_curve=3000, + penalty_upper_rule_curve=3000, + penalty_final_level=3000, + force_final_level=True, + ) + + all_reservoirs = MultiStockManagement( + [reservoir_management_1, reservoir_management_2] + ) + + x_1 = np.linspace(0, reservoir_1.capacity, num=5) + x_2 = np.linspace(0, reservoir_2.capacity, num=5) + + try: + + vb, lb, ub = calculate_bellman_value_directly( + param=param, + multi_stock_management=all_reservoirs, + output_path="test_data/two_nodes", + X={"area_1": x_1, "area_2": x_2}, + univariate=True, + ) + assert False + except AssertionError: + assert True diff --git a/tests/test_bellman_value_precalculated_reward_two_stocks.py b/tests/test_bellman_value_precalculated_reward_two_stocks.py index f3e4841..c27c5ff 100644 --- a/tests/test_bellman_value_precalculated_reward_two_stocks.py +++ b/tests/test_bellman_value_precalculated_reward_two_stocks.py @@ -1,10 +1,10 @@ import numpy as np import pytest -from calculate_reward_and_bellman_values import MultiStockManagement from functions_iterative import ReservoirManagement, TimeScenarioParameter from multi_stock_bellman_value_calculation import * from read_antares_data import Reservoir +from reservoir_management import MultiStockManagement def test_bellman_value_precalculated_multi_stock() -> None: diff --git a/tests/test_fast_uv_two_stocks.py b/tests/test_fast_uv_two_stocks.py index 5ee3eb9..6a290e7 100644 --- a/tests/test_fast_uv_two_stocks.py +++ b/tests/test_fast_uv_two_stocks.py @@ -1,10 +1,10 @@ import numpy as np import pytest -from calculate_reward_and_bellman_values import MultiStockManagement from functions_iterative import ReservoirManagement, TimeScenarioParameter from multi_stock_bellman_value_calculation import generate_fast_uvs_v2 from read_antares_data import Reservoir +from reservoir_management import MultiStockManagement def test_fast_usage_values_multi_stock() -> None: diff --git a/tests/test_iterative_two_stocks.py b/tests/test_iterative_two_stocks.py index ca27181..1970f22 100644 --- a/tests/test_iterative_two_stocks.py +++ b/tests/test_iterative_two_stocks.py @@ -1,10 +1,10 @@ import numpy as np import pytest -from calculate_reward_and_bellman_values import MultiStockManagement from functions_iterative import ReservoirManagement, TimeScenarioParameter from multi_stock_bellman_value_calculation import * from read_antares_data import Reservoir +from reservoir_management import MultiStockManagement def test_bellman_value_iterative_method() -> None: diff --git a/tests/test_optimization_problem_two_stock.py b/tests/test_optimization_problem_two_stock.py index cfa2060..8666a90 100644 --- a/tests/test_optimization_problem_two_stock.py +++ b/tests/test_optimization_problem_two_stock.py @@ -1,10 +1,10 @@ import ortools.linear_solver.pywraplp as pywraplp import pytest -from calculate_reward_and_bellman_values import MultiStockManagement from functions_iterative import ReservoirManagement, TimeScenarioParameter from optimization import AntaresProblem, Basis from read_antares_data import Reservoir +from reservoir_management import MultiStockManagement def test_create_weekly_problem_with_two_stocks() -> None: diff --git a/tests/test_read_optimization_problem.py b/tests/test_read_optimization_problem.py index 9ea30fb..4976d67 100644 --- a/tests/test_read_optimization_problem.py +++ b/tests/test_read_optimization_problem.py @@ -5,7 +5,6 @@ from calculate_reward_and_bellman_values import ( BellmanValueCalculation, MultiStockBellmanValueCalculation, - MultiStockManagement, RewardApproximation, ) from estimation import PieceWiseLinearInterpolator, UniVariateEstimator @@ -16,6 +15,7 @@ ) from optimization import AntaresProblem, Basis from read_antares_data import Reservoir +from reservoir_management import MultiStockManagement def test_create_and_modify_weekly_problem() -> None: diff --git a/tests/test_upper_bound.py b/tests/test_upper_bound.py index 3011828..a84bfa2 100644 --- a/tests/test_upper_bound.py +++ b/tests/test_upper_bound.py @@ -2,10 +2,7 @@ import ortools.linear_solver.pywraplp as pywraplp import pytest -from calculate_reward_and_bellman_values import ( - MultiStockBellmanValueCalculation, - MultiStockManagement, -) +from calculate_reward_and_bellman_values import MultiStockBellmanValueCalculation from estimation import PieceWiseLinearInterpolator, UniVariateEstimator from functions_iterative import ( BellmanValueCalculation, @@ -17,6 +14,7 @@ ) from optimization import AntaresProblem from read_antares_data import Reservoir +from reservoir_management import MultiStockManagement bellman_values = np.array( [