diff --git a/src/andromede/simulation/__init__.py b/src/andromede/simulation/__init__.py index afdb1dc4..23e0a855 100644 --- a/src/andromede/simulation/__init__.py +++ b/src/andromede/simulation/__init__.py @@ -15,6 +15,7 @@ build_benders_decomposed_problem, ) from .optimization import BlockBorderManagement, OptimizationProblem, build_problem -from .output_values import OutputValues +from .output_values import BendersSolution, OutputValues +from .runner import BendersRunner, MergeMPSRunner from .strategy import MergedProblemStrategy, ModelSelectionStrategy from .time_block import TimeBlock diff --git a/src/andromede/simulation/benders_decomposed.py b/src/andromede/simulation/benders_decomposed.py index e9dd9a6a..08d718f1 100644 --- a/src/andromede/simulation/benders_decomposed.py +++ b/src/andromede/simulation/benders_decomposed.py @@ -15,18 +15,20 @@ with Benders solver related functions """ -import json -import os import pathlib -import subprocess -import sys -from typing import Any, Dict, List +from typing import Any, Dict, List, Optional from andromede.simulation.optimization import ( BlockBorderManagement, OptimizationProblem, build_problem, ) +from andromede.simulation.output_values import ( + BendersDecomposedSolution, + BendersMergedSolution, + BendersSolution, +) +from andromede.simulation.runner import BendersRunner, MergeMPSRunner from andromede.simulation.strategy import ( InvestmentProblemStrategy, OperationalProblemStrategy, @@ -34,7 +36,7 @@ from andromede.simulation.time_block import TimeBlock from andromede.study.data import DataBase from andromede.study.network import Network -from andromede.utils import serialize +from andromede.utils import read_json, serialize, serialize_json class BendersDecomposedProblem: @@ -45,12 +47,28 @@ class BendersDecomposedProblem: master: OptimizationProblem subproblems: List[OptimizationProblem] + emplacement: pathlib.Path + output_path: pathlib.Path + + solution: Optional[BendersSolution] + is_merged: bool + def __init__( - self, master: OptimizationProblem, subproblems: List[OptimizationProblem] + self, + master: OptimizationProblem, + subproblems: List[OptimizationProblem], + emplacement: str = "outputs/lp", + output_path: str = "expansion", ) -> None: self.master = master self.subproblems = subproblems + self.emplacement = pathlib.Path(emplacement) + self.output_path = pathlib.Path(output_path) + + self.solution = None + self.is_merged = False + def export_structure(self) -> str: """ Write the structure.txt file @@ -111,8 +129,8 @@ def export_options( "BOUND_ALPHA": True, "SEPARATION_PARAM": 0.5, "BATCH_SIZE": 0, - "JSON_FILE": "output/xpansion/out.json", - "LAST_ITERATION_JSON_FILE": "output/xpansion/last_iteration.json", + "JSON_FILE": f"{self.output_path}/out.json", + "LAST_ITERATION_JSON_FILE": f"{self.output_path}/last_iteration.json", "MASTER_FORMULATION": "integer", "SOLVER_NAME": solver_name, "TIME_LIMIT": 1_000_000_000_000, @@ -125,62 +143,70 @@ def export_options( def prepare( self, *, - path: str = "outputs/lp", solver_name: str = "XPRESS", log_level: int = 0, is_debug: bool = False, ) -> None: - directory = pathlib.Path(path) - serialize("master.mps", self.master.export_as_mps(), directory) - serialize("subproblem.mps", self.subproblems[0].export_as_mps(), directory) - serialize("structure.txt", self.export_structure(), directory) - serialize( + serialize("master.mps", self.master.export_as_mps(), self.emplacement) + for subproblem in self.subproblems: + serialize( + f"{subproblem.name}.mps", subproblem.export_as_mps(), self.emplacement + ) + serialize("structure.txt", self.export_structure(), self.emplacement) + serialize_json( "options.json", - json.dumps( - self.export_options(solver_name=solver_name, log_level=log_level), - indent=4, - ), - directory, + self.export_options(solver_name=solver_name, log_level=log_level), + self.emplacement, ) if is_debug: - serialize("master.lp", self.master.export_as_lp(), directory) - serialize("subproblem.lp", self.subproblems[0].export_as_lp(), directory) + serialize("master.lp", self.master.export_as_lp(), self.emplacement) + for subproblem in self.subproblems: + serialize( + f"{subproblem.name}.lp", subproblem.export_as_lp(), self.emplacement + ) + + def read_solution(self) -> None: + try: + data = read_json("out.json", self.emplacement / self.output_path) + + except FileNotFoundError: + # TODO For now, it will return as if nothing is wrong + # modify it with runner's run + print("Return without reading it for now") + return + + if self.is_merged: + self.solution = BendersMergedSolution(data) + else: + self.solution = BendersDecomposedSolution(data) def run( self, *, - path: str = "outputs/lp", solver_name: str = "XPRESS", log_level: int = 0, + should_merge: bool = False, ) -> bool: - self.prepare(path=path, solver_name=solver_name, log_level=log_level) - root_dir = pathlib.Path().cwd() - path_to_benders = root_dir / "bin" / "benders" - - if not path_to_benders.is_file(): - # TODO Maybe a more robust check and/or return value? - # For now, it won't look anywhere else because a new - # architecture should be discussed - print(f"{path_to_benders} executable not found. Returning True") - return True + self.prepare(solver_name=solver_name, log_level=log_level) - os.chdir(path) - res = subprocess.run( - [path_to_benders, "options.json"], - stdout=sys.stdout, - stderr=subprocess.DEVNULL, # TODO For now, to avoid the "Invalid MIT-MAGIC-COOKIE-1 key" error - shell=False, - ) - os.chdir(root_dir) + if not should_merge: + return_code = BendersRunner(self.emplacement).run() + else: + self.is_merged = True + return_code = MergeMPSRunner(self.emplacement).run() - return res.returncode == 0 + if return_code == 0: + self.read_solution() + return True + else: + return False def build_benders_decomposed_problem( network: Network, database: DataBase, - block: TimeBlock, + blocks: List[TimeBlock], scenarios: int, *, border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, @@ -196,8 +222,10 @@ def build_benders_decomposed_problem( master = build_problem( network, database, - block, - scenarios, + null_time_block := TimeBlock( # Not necessary for master, but list must be non-empty + 0, [0] + ), + null_scenario := 0, # Not necessary for master problem_name="master", border_management=border_management, solver_id=solver_id, @@ -205,15 +233,19 @@ def build_benders_decomposed_problem( ) # Benders Decomposed Sub-problems - subproblem = build_problem( - network, - database, - block, - scenarios, - problem_name="subproblem", - border_management=border_management, - solver_id=solver_id, - problem_strategy=OperationalProblemStrategy(), - ) + subproblems = [] + for block in blocks: + subproblems.append( + build_problem( + network, + database, + block, + scenarios, + problem_name=f"subproblem_{block.id}", + border_management=border_management, + solver_id=solver_id, + problem_strategy=OperationalProblemStrategy(), + ) + ) - return BendersDecomposedProblem(master, [subproblem]) + return BendersDecomposedProblem(master, subproblems) diff --git a/src/andromede/simulation/optimization.py b/src/andromede/simulation/optimization.py index 7e26dd38..35c50e4a 100644 --- a/src/andromede/simulation/optimization.py +++ b/src/andromede/simulation/optimization.py @@ -631,7 +631,7 @@ def make_constraint( Adds constraint to the solver. """ solver_constraints = {} - constraint_name = data.name + constraint_name = f"{data.name}_t{block_timestep}_s{scenario}" for instance in range(instances): if instances > 1: constraint_name += f"_{instance}" @@ -753,7 +753,7 @@ def _create_variables(self) -> None: solver_var = self.solver.NumVar( lower_bound, upper_bound, - f"{component.id}_{model_var.name}_{block_timestep}_{scenario}", + f"{component.id}_{model_var.name}_t{block_timestep}_s{scenario}", ) component_context.add_variable( block_timestep, scenario, model_var.name, solver_var diff --git a/src/andromede/simulation/output_values.py b/src/andromede/simulation/output_values.py index d875bb8f..cf94e7c0 100644 --- a/src/andromede/simulation/output_values.py +++ b/src/andromede/simulation/output_values.py @@ -15,7 +15,7 @@ """ import math from dataclasses import dataclass, field -from typing import Dict, List, Mapping, Optional, Tuple, TypeVar, Union, cast +from typing import Any, Dict, List, Mapping, Optional, Tuple, TypeVar, Union, cast from andromede.simulation.optimization import OptimizationProblem from andromede.study.data import TimeScenarioIndex @@ -247,3 +247,107 @@ def _are_mappings_close( ) else: return True + + +@dataclass(frozen=True) +class BendersSolution: + data: Dict[str, Any] + + def __eq__(self, other: object) -> bool: + if not isinstance(other, BendersSolution): + return NotImplemented + return ( + self.overall_cost == other.overall_cost + and self.candidates == other.candidates + ) + + def is_close( + self, + other: "BendersSolution", + *, + rel_tol: float = 1.0e-9, + abs_tol: float = 0.0, + ) -> bool: + return ( + math.isclose( + self.overall_cost, other.overall_cost, abs_tol=abs_tol, rel_tol=rel_tol + ) + and self.candidates.keys() == other.candidates.keys() + and all( + math.isclose( + self.candidates[key], + other.candidates[key], + rel_tol=rel_tol, + abs_tol=abs_tol, + ) + for key in self.candidates + ) + ) + + def __str__(self) -> str: + lpad = 30 + rpad = 12 + + string = "Benders' solution:\n" + string += f"{'Overall cost':<{lpad}} : {self.overall_cost:>{rpad}}\n" + string += f"{'Investment cost':<{lpad}} : {self.investment_cost:>{rpad}}\n" + string += f"{'Operational cost':<{lpad}} : {self.operational_cost:>{rpad}}\n" + string += "-" * (lpad + rpad + 3) + "\n" + for candidate, investment in self.candidates.items(): + string += f"{candidate:<{lpad}} : {investment:>{rpad}}\n" + + return string + + @property + def investment_cost(self) -> float: + return self.data["solution"]["investment_cost"] + + @property + def operational_cost(self) -> float: + return self.data["solution"]["operational_cost"] + + @property + def overall_cost(self) -> float: + return self.data["solution"]["overall_cost"] + + @property + def candidates(self) -> Dict[str, float]: + return self.data["solution"]["values"] + + @property + def status(self) -> str: + return self.data["solution"]["problem_status"] + + @property + def absolute_gap(self) -> float: + return self.data["solution"]["optimality_gap"] + + @property + def relative_gap(self) -> float: + return self.data["solution"]["relative_gap"] + + @property + def stopping_criterion(self) -> str: + return self.data["solution"]["stopping_criterion"] + + +@dataclass(frozen=True, eq=False) +class BendersMergedSolution(BendersSolution): + @property + def lower_bound(self) -> float: + return self.data["solution"]["lb"] + + @property + def upper_bound(self) -> float: + return self.data["solution"]["ub"] + + +@dataclass(frozen=True, eq=False) +class BendersDecomposedSolution(BendersSolution): + @property + def nb_iterations(self) -> int: + return self.data["solution"]["iteration"] + + @property + def duration(self) -> float: + return self.data["run_duration"] diff --git a/src/andromede/simulation/runner.py b/src/andromede/simulation/runner.py new file mode 100644 index 00000000..ac52a04f --- /dev/null +++ b/src/andromede/simulation/runner.py @@ -0,0 +1,66 @@ +# 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 os +import pathlib +import subprocess +import sys +from typing import List + + +class CommandRunner: + def __init__( + self, + binary_path: pathlib.Path, + list_arguments: List[str], + emplacement: pathlib.Path, + ) -> None: + self.current_dir: pathlib.Path = pathlib.Path().cwd() + self.command: pathlib.Path = binary_path + self.emplacement: pathlib.Path = emplacement + self.arguments: List[str] = list_arguments + + def check_command(self) -> bool: + if not self.command.is_file(): + print(f"{self.command} executable not found") + return False + return True + + def run(self) -> int: + if not self.check_command(): + # TODO For now, it will return 0 as if nothing is wrong + # eventually if should return an error + # maybe wait when we separate unit tests from integration tests + # modify with bender_decomposed's read_solution + print("Return code 0 for now") + return 0 + + os.chdir(self.emplacement) + res = subprocess.run( + [self.current_dir / self.command, *self.arguments], + stdout=sys.stdout, + stderr=subprocess.DEVNULL, # TODO For now, to avoid the "Invalid MIT-MAGIC-COOKIE-1 key" error + shell=False, + ) + os.chdir(self.current_dir) + + return res.returncode + + +class BendersRunner(CommandRunner): + def __init__(self, emplacement: pathlib.Path) -> None: + super().__init__(pathlib.Path("bin/benders"), ["options.json"], emplacement) + + +class MergeMPSRunner(CommandRunner): + def __init__(self, emplacement: pathlib.Path) -> None: + super().__init__(pathlib.Path("bin/merge_mps"), ["options.json"], emplacement) diff --git a/src/andromede/utils.py b/src/andromede/utils.py index 1b5c00b8..baa32ff5 100644 --- a/src/andromede/utils.py +++ b/src/andromede/utils.py @@ -13,8 +13,9 @@ """ Module for technical utilities. """ +import json import pathlib -from typing import Any, Callable, Dict, Generic, Optional, TypeVar +from typing import Any, Callable, Dict, Optional, TypeVar T = TypeVar("T") K = TypeVar("K") @@ -55,3 +56,15 @@ def serialize(filename: str, message: str, path: pathlib.Path) -> None: with file: file.write(message) + + +def serialize_json( + filename: str, message: Dict[str, Any], path: pathlib.Path, indentation: int = 4 +) -> None: + serialize(filename, json.dumps(message, indent=indentation), path) + + +def read_json(filename: str, path: pathlib.Path) -> Dict[str, Any]: + with (path / filename).open() as file: + data = json.load(file) + return data diff --git a/tests/integration/test_benders_decomposed.py b/tests/integration/test_benders_decomposed.py index b8e835ec..38cc06b0 100644 --- a/tests/integration/test_benders_decomposed.py +++ b/tests/integration/test_benders_decomposed.py @@ -33,7 +33,11 @@ model, ) from andromede.model.model import PortFieldDefinition, PortFieldId -from andromede.simulation import TimeBlock, build_benders_decomposed_problem +from andromede.simulation import ( + BendersSolution, + TimeBlock, + build_benders_decomposed_problem, +) from andromede.study import ( Component, ConstantData, @@ -41,6 +45,8 @@ Network, Node, PortRef, + TimeIndex, + TimeSeriesData, create_component, ) @@ -161,14 +167,33 @@ def cluster_candidate(discrete_candidate: Model) -> Component: return cluster -def test_benders_decomposed_single_time_step_single_scenario( +def test_benders_decomposed_integration( generator: Component, candidate: Component, cluster_candidate: Component, ) -> None: """ - Same test as before but this time we separate master/subproblem and - export the problems in MPS format to be solved by the Benders solver in Xpansion + Simple generation expansion problem on one node, one timestep and one scenario + but this time with two candidates: one thermal cluster and one wind cluster. + We separate master/subproblem and export the problems in MPS format to be solved by the Benders solver in Xpansion + + Demand = 400 + Generator : P_max : 200, Cost : 45 + Unsupplied energy : Cost : 501 + + -> 200 of unsupplied energy + -> Total cost without investment = 45 * 200 + 501 * 200 = 109_200 + + Single candidate : Invest cost : 490 / MW; Prod cost : 10 + Cluster candidate : Invest cost : 200 / MW; Prod cost : 10; Nb of discrete thresholds: 10; Prod per threshold: 10 + + Optimal investment : 100 MW (Cluster) + 100 MW (Single) + + -> Optimal cost = 490 * 100 + 10 * 100 (Single) + + 200 * 100 + 10 * 100 (Cluster) + + 45 * 200 (Generator) + = 69_000 + 11_000 + = 80_000 """ database = DataBase() @@ -205,6 +230,113 @@ def test_benders_decomposed_single_time_step_single_scenario( scenarios = 1 xpansion = build_benders_decomposed_problem( - network, database, TimeBlock(1, [0]), scenarios + network, database, [TimeBlock(1, [0])], scenarios ) + + data = { + "solution": { + "overall_cost": 80_000, + "values": {"CAND_p_max_t0_s0": 100, "DISCRETE_p_max_t0_s0": 100}, + } + } + solution = BendersSolution(data) + + assert xpansion.run() + decomposed_solution = xpansion.solution + if decomposed_solution is not None: # For mypy only + assert decomposed_solution.is_close( + solution + ), f"Solution differs from expected: {decomposed_solution}" + + assert xpansion.run(should_merge=True) + merged_solution = xpansion.solution + if merged_solution is not None: # For mypy only + assert merged_solution.is_close( + solution + ), f"Solution differs from expected: {merged_solution}" + + +def test_benders_decomposed_multi_time_block_single_scenario( + generator: Component, + candidate: Component, +) -> None: + """ + Simple generation xpansion problem on one node. Two time blocks with one timestep each, + one scenario, one thermal cluster candidate. + + Demand = [200, 300] + Generator : P_max : 200, Cost : 40 + Unsupplied energy : Cost : 501 + + -> [0, 100] of unsupplied energy + -> Total cost without investment = (200 * 40) + (200 * 40 + 100 * 501) = 66_100 + + Candidate : Invest cost : 480 / MW, Prod cost : 10 + + Optimal investment : 100 MW + + -> Optimal cost = 480 * 100 (investment) + + 10 * 100 + 40 * 100 (operational - time block 1) + + 10 * 100 + 40 * 200 (operational - time block 2) + = 62_000 + + """ + + data = {} + data[TimeIndex(0)] = 200 + data[TimeIndex(1)] = 300 + + demand_data = TimeSeriesData(time_series=data) + + database = DataBase() + database.add_data("D", "demand", demand_data) + + database.add_data("N", "spillage_cost", ConstantData(1)) + database.add_data("N", "ens_cost", ConstantData(501)) + + database.add_data("G1", "p_max", ConstantData(200)) + database.add_data("G1", "cost", ConstantData(40)) + + database.add_data("CAND", "op_cost", ConstantData(10)) + database.add_data("CAND", "invest_cost", ConstantData(480)) + + demand = create_component( + model=DEMAND_MODEL, + id="D", + ) + + node = Node(model=NODE_WITH_SPILL_AND_ENS_MODEL, id="N") + network = Network("test") + network.add_node(node) + network.add_component(demand) + network.add_component(generator) + network.add_component(candidate) + network.connect(PortRef(demand, "balance_port"), PortRef(node, "balance_port")) + network.connect(PortRef(generator, "balance_port"), PortRef(node, "balance_port")) + network.connect(PortRef(candidate, "balance_port"), PortRef(node, "balance_port")) + + scenarios = 1 + + xpansion = build_benders_decomposed_problem( + network, + database, + [TimeBlock(1, [0]), TimeBlock(2, [1])], + scenarios, + ) + + data = { + "solution": { + "overall_cost": 62_000, + "values": { + "CAND_p_max_t0_s0": 100, + }, + } + } + solution = BendersSolution(data) + assert xpansion.run() + decomposed_solution = xpansion.solution + if decomposed_solution is not None: # For mypy only + assert decomposed_solution.is_close( + solution + ), f"Solution differs from expected: {decomposed_solution}"