From 9fc446479def5d6f31a261281b59e8df4cf3f831 Mon Sep 17 00:00:00 2001 From: ianmnz <99749877+ianmnz@users.noreply.github.com> Date: Mon, 26 Feb 2024 15:48:16 +0100 Subject: [PATCH] Benders decomposition for investment problems (#12) - Add benders-decomposed problem resolution, which will call benders tool from antares-xpansion. - Variables, constraints and objectives can be tagged as investment or operational. - Problem building from the data can now be customized with a strategy to select which variables etc. to include. This strategy has 2 implementations: operational and investment. --- .gitignore | 2 + src/andromede/libs/standard.py | 58 ++- src/andromede/model/__init__.py | 5 +- src/andromede/model/common.py | 28 ++ src/andromede/model/constraint.py | 5 + src/andromede/model/model.py | 44 +- src/andromede/model/parameter.py | 14 +- src/andromede/model/variable.py | 41 +- src/andromede/simulation/__init__.py | 7 +- .../simulation/benders_decomposed.py | 219 ++++++++++ src/andromede/simulation/optimization.py | 382 ++++++++++-------- src/andromede/simulation/output_values.py | 4 +- src/andromede/simulation/strategy.py | 82 ++++ src/andromede/utils.py | 14 + tests/andromede/test_andromede.py | 21 +- tests/andromede/test_data.py | 8 +- tests/andromede/test_electrolyzer.py | 4 +- tests/andromede/test_model.py | 2 +- tests/andromede/test_output_values.py | 10 +- tests/andromede/test_performance.py | 6 +- tests/andromede/test_xpansion.py | 230 ++++++++++- 21 files changed, 928 insertions(+), 258 deletions(-) create mode 100644 src/andromede/model/common.py create mode 100644 src/andromede/simulation/benders_decomposed.py create mode 100644 src/andromede/simulation/strategy.py diff --git a/.gitignore b/.gitignore index cb215344..e1b8586c 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,5 @@ venv .env .coverage coverage.xml +bin +outputs diff --git a/src/andromede/libs/standard.py b/src/andromede/libs/standard.py index 856cbcb5..08135f56 100644 --- a/src/andromede/libs/standard.py +++ b/src/andromede/libs/standard.py @@ -41,6 +41,28 @@ ], ) +NODE_WITH_SPILL_AND_ENS_MODEL = model( + id="NODE_WITH_SPILL_AND_ENS_MODEL", + parameters=[float_parameter("spillage_cost"), float_parameter("ens_cost")], + variables=[ + float_variable("spillage", lower_bound=literal(0)), + float_variable("unsupplied_energy", lower_bound=literal(0)), + ], + ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")], + binding_constraints=[ + Constraint( + name="Balance", + expression=port_field("balance_port", "flow").sum_connections() + == var("spillage") - var("unsupplied_energy"), + ) + ], + objective_operational_contribution=( + param("spillage_cost") * var("spillage") + + param("ens_cost") * var("unsupplied_energy") + ) + .sum() + .expec(), +) """ A standard model for a linear cost generation, limited by a maximum generation. """ @@ -63,7 +85,9 @@ name="Max generation", expression=var("generation") <= param("p_max") ), ], - objective_contribution=(param("cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("cost") * var("generation")) + .sum() + .expec(), ) """ @@ -133,7 +157,9 @@ lower_bound=literal(0), ), # To test both ways of setting constraints ], - objective_contribution=(param("cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("cost") * var("generation")) + .sum() + .expec(), ) # For now, no starting cost @@ -159,17 +185,17 @@ "nb_on", lower_bound=literal(0), upper_bound=param("nb_units_max"), - structural_type=ANTICIPATIVE_TIME_VARYING, + structure=ANTICIPATIVE_TIME_VARYING, ), int_variable( "nb_stop", lower_bound=literal(0), - structural_type=ANTICIPATIVE_TIME_VARYING, + structure=ANTICIPATIVE_TIME_VARYING, ), int_variable( "nb_start", lower_bound=literal(0), - structural_type=ANTICIPATIVE_TIME_VARYING, + structure=ANTICIPATIVE_TIME_VARYING, ), ], ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")], @@ -208,7 +234,9 @@ ) # It also works by writing ExpressionRange(-param("d_min_down") + 1, 0) as ExpressionRange's __post_init__ wraps integers to literal nodes. However, MyPy does not seem to infer that ExpressionRange's attributes are necessarily of ExpressionNode type and raises an error if the arguments in the constructor are integer (whereas it runs correctly), this why we specify it here with literal(0) instead of 0. ], - objective_contribution=(param("cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("cost") * var("generation")) + .sum() + .expec(), ) # Same model as previous one, except that starting/stopping variables are now non anticipative @@ -234,17 +262,17 @@ "nb_on", lower_bound=literal(0), upper_bound=param("nb_units_max"), - structural_type=NON_ANTICIPATIVE_TIME_VARYING, + structure=NON_ANTICIPATIVE_TIME_VARYING, ), int_variable( "nb_stop", lower_bound=literal(0), - structural_type=NON_ANTICIPATIVE_TIME_VARYING, + structure=NON_ANTICIPATIVE_TIME_VARYING, ), int_variable( "nb_start", lower_bound=literal(0), - structural_type=NON_ANTICIPATIVE_TIME_VARYING, + structure=NON_ANTICIPATIVE_TIME_VARYING, ), ], ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")], @@ -282,7 +310,9 @@ <= param("nb_units_max").shift(-param("d_min_down")) - var("nb_on"), ), ], - objective_contribution=(param("cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("cost") * var("generation")) + .sum() + .expec(), ) SPILLAGE_MODEL = model( @@ -296,7 +326,7 @@ definition=-var("spillage"), ) ], - objective_contribution=(param("cost") * var("spillage")).sum().expec(), + objective_operational_contribution=(param("cost") * var("spillage")).sum().expec(), ) UNSUPPLIED_ENERGY_MODEL = model( @@ -310,7 +340,9 @@ definition=var("unsupplied_energy"), ) ], - objective_contribution=(param("cost") * var("unsupplied_energy")).sum().expec(), + objective_operational_contribution=(param("cost") * var("unsupplied_energy")) + .sum() + .expec(), ) # Simplified model @@ -357,5 +389,5 @@ == param("inflows"), ), ], - objective_contribution=literal(0), # Implcitement nul ? + objective_operational_contribution=literal(0), # Implcitement nul ? ) diff --git a/src/andromede/model/__init__.py b/src/andromede/model/__init__.py index db22b519..6aee655b 100644 --- a/src/andromede/model/__init__.py +++ b/src/andromede/model/__init__.py @@ -10,8 +10,9 @@ # # This file is part of the Antares project. +from .common import ProblemContext, ValueType from .constraint import Constraint from .model import Model, ModelPort, model -from .parameter import Parameter, ParameterValueType, float_parameter, int_parameter +from .parameter import Parameter, float_parameter, int_parameter from .port import PortField, PortType -from .variable import Variable, VariableValueType, float_variable, int_variable +from .variable import Variable, float_variable, int_variable diff --git a/src/andromede/model/common.py b/src/andromede/model/common.py new file mode 100644 index 00000000..07ebdcb4 --- /dev/null +++ b/src/andromede/model/common.py @@ -0,0 +1,28 @@ +# 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. + +""" +Module for common classes used in models. +""" +from enum import Enum + + +class ValueType(Enum): + FLOAT = "FLOAT" + INTEGER = "INTEGER" + # Needs more ? + + +class ProblemContext(Enum): + OPERATIONAL = 0 + INVESTMENT = 1 + COUPLING = 2 diff --git a/src/andromede/model/constraint.py b/src/andromede/model/constraint.py index 2418f595..9e7c6050 100644 --- a/src/andromede/model/constraint.py +++ b/src/andromede/model/constraint.py @@ -20,6 +20,7 @@ literal, ) from andromede.expression.print import print_expr +from andromede.model.common import ProblemContext class Constraint: @@ -33,6 +34,7 @@ class Constraint: expression: ExpressionNode lower_bound: ExpressionNode upper_bound: ExpressionNode + context: ProblemContext def __init__( self, @@ -40,8 +42,11 @@ def __init__( expression: ExpressionNode, lower_bound: Optional[ExpressionNode] = None, upper_bound: Optional[ExpressionNode] = None, + context: ProblemContext = ProblemContext.OPERATIONAL, ) -> None: self.name = name + self.context = context + if isinstance(expression, ComparisonNode): if lower_bound is not None or upper_bound is not None: raise ValueError( diff --git a/src/andromede/model/model.py b/src/andromede/model/model.py index ba54ac31..3997f43d 100644 --- a/src/andromede/model/model.py +++ b/src/andromede/model/model.py @@ -78,6 +78,23 @@ def get_component_variable_structure( return Provider() +def _is_objective_contribution_valid( + model: "Model", objective_contribution: ExpressionNode +) -> bool: + if not is_linear(objective_contribution): + raise ValueError("Objective contribution must be a linear expression.") + + data_structure_provider = _make_structure_provider(model) + objective_structure = compute_indexation( + objective_contribution, data_structure_provider + ) + + if objective_structure != IndexingStructure(time=False, scenario=False): + raise ValueError("Objective contribution should be a real-valued expression.") + # TODO: We should also check that the number of instances is equal to 1, but this would require a linearization here, do not want to do that for now... + return True + + @dataclass(frozen=True) class ModelPort: """ @@ -129,26 +146,23 @@ class Model: inter_block_dyn: bool = False parameters: Dict[str, Parameter] = field(default_factory=dict) variables: Dict[str, Variable] = field(default_factory=dict) - objective_contribution: Optional[ExpressionNode] = None + objective_operational_contribution: Optional[ExpressionNode] = None + objective_investment_contribution: Optional[ExpressionNode] = None ports: Dict[str, ModelPort] = field(default_factory=dict) # key = port name port_fields_definitions: Dict[PortFieldId, PortFieldDefinition] = field( default_factory=dict ) def __post_init__(self) -> None: - if self.objective_contribution: - if not is_linear(self.objective_contribution): - raise ValueError("Objective contribution must be a linear expression.") + if self.objective_operational_contribution: + _is_objective_contribution_valid( + self, self.objective_operational_contribution + ) - data_structure_provider = _make_structure_provider(self) - objective_structure = compute_indexation( - self.objective_contribution, data_structure_provider + if self.objective_investment_contribution: + _is_objective_contribution_valid( + self, self.objective_investment_contribution ) - if objective_structure != IndexingStructure(time=False, scenario=False): - raise ValueError( - "Objective contribution should be a real-valued expression." - ) - # TODO: We should also check that the number of instances is equal to 1, but this would require a linearization here, do not want to do that for now... for definition in self.port_fields_definitions.values(): port_name = definition.port_field.port_name @@ -176,7 +190,8 @@ def model( binding_constraints: Optional[Iterable[Constraint]] = None, parameters: Optional[Iterable[Parameter]] = None, variables: Optional[Iterable[Variable]] = None, - objective_contribution: Optional[ExpressionNode] = None, + objective_operational_contribution: Optional[ExpressionNode] = None, + objective_investment_contribution: Optional[ExpressionNode] = None, inter_block_dyn: bool = False, ports: Optional[Iterable[ModelPort]] = None, port_fields_definitions: Optional[Iterable[PortFieldDefinition]] = None, @@ -202,7 +217,8 @@ def model( else {}, parameters={p.name: p for p in parameters} if parameters else {}, variables={v.name: v for v in variables} if variables else {}, - objective_contribution=objective_contribution, + objective_operational_contribution=objective_operational_contribution, + objective_investment_contribution=objective_investment_contribution, inter_block_dyn=inter_block_dyn, ports=existing_port_names, port_fields_definitions={d.port_field: d for d in port_fields_definitions} diff --git a/src/andromede/model/parameter.py b/src/andromede/model/parameter.py index 2ec2dbb0..745d5a40 100644 --- a/src/andromede/model/parameter.py +++ b/src/andromede/model/parameter.py @@ -11,15 +11,9 @@ # This file is part of the Antares project. from dataclasses import dataclass -from enum import Enum from andromede.expression.indexing_structure import IndexingStructure - - -class ParameterValueType(Enum): - FLOAT = "FLOAT" - INTEGER = "INTEGER" - # Needs more ? +from andromede.model.common import ValueType @dataclass(frozen=True) @@ -31,7 +25,7 @@ class Parameter: """ name: str - type: ParameterValueType + type: ValueType structure: IndexingStructure @@ -39,11 +33,11 @@ def int_parameter( name: str, structure: IndexingStructure = IndexingStructure(True, True), ) -> Parameter: - return Parameter(name, ParameterValueType.INTEGER, structure) + return Parameter(name, ValueType.INTEGER, structure) def float_parameter( name: str, structure: IndexingStructure = IndexingStructure(True, True), ) -> Parameter: - return Parameter(name, ParameterValueType.FLOAT, structure) + return Parameter(name, ValueType.FLOAT, structure) diff --git a/src/andromede/model/variable.py b/src/andromede/model/variable.py index 5ddc2d93..7f65b2c4 100644 --- a/src/andromede/model/variable.py +++ b/src/andromede/model/variable.py @@ -10,19 +10,13 @@ # # This file is part of the Antares project. -from dataclasses import dataclass, field -from enum import Enum +from dataclasses import dataclass from typing import Optional from andromede.expression import ExpressionNode from andromede.expression.degree import is_constant from andromede.expression.indexing_structure import IndexingStructure - - -class VariableValueType(Enum): - FLOAT = "FLOAT" - INTEGER = "INTEGER" - # Needs more ? +from andromede.model.common import ProblemContext, ValueType @dataclass @@ -32,43 +26,36 @@ class Variable: """ name: str - data_type: VariableValueType + data_type: ValueType lower_bound: Optional[ExpressionNode] upper_bound: Optional[ExpressionNode] - structure: IndexingStructure = field(default=IndexingStructure(True, True)) + structure: IndexingStructure + context: ProblemContext def __post_init__(self) -> None: if self.lower_bound and not is_constant(self.lower_bound): raise ValueError("Lower bounds of variables must be constant") if self.upper_bound and not is_constant(self.upper_bound): - raise ValueError("Lower bounds of variables must be constant") + raise ValueError("Upper bounds of variables must be constant") def int_variable( name: str, lower_bound: Optional[ExpressionNode] = None, upper_bound: Optional[ExpressionNode] = None, - structural_type: Optional[IndexingStructure] = None, + structure: IndexingStructure = IndexingStructure(True, True), + context: ProblemContext = ProblemContext.OPERATIONAL, ) -> Variable: - # Dirty if/else just for MyPy - if structural_type is None: - return Variable(name, VariableValueType.INTEGER, lower_bound, upper_bound) - else: - return Variable( - name, VariableValueType.INTEGER, lower_bound, upper_bound, structural_type - ) + return Variable( + name, ValueType.INTEGER, lower_bound, upper_bound, structure, context + ) def float_variable( name: str, lower_bound: Optional[ExpressionNode] = None, upper_bound: Optional[ExpressionNode] = None, - structure: Optional[IndexingStructure] = None, + structure: IndexingStructure = IndexingStructure(True, True), + context: ProblemContext = ProblemContext.OPERATIONAL, ) -> Variable: - # Dirty if/else just for MyPy - if structure is None: - return Variable(name, VariableValueType.FLOAT, lower_bound, upper_bound) - else: - return Variable( - name, VariableValueType.FLOAT, lower_bound, upper_bound, structure - ) + return Variable(name, ValueType.FLOAT, lower_bound, upper_bound, structure, context) diff --git a/src/andromede/simulation/__init__.py b/src/andromede/simulation/__init__.py index b7c99b3b..afdb1dc4 100644 --- a/src/andromede/simulation/__init__.py +++ b/src/andromede/simulation/__init__.py @@ -10,6 +10,11 @@ # # This file is part of the Antares project. -from .optimization import BlockBorderManagement, SolverAndContext, build_problem +from .benders_decomposed import ( + BendersDecomposedProblem, + build_benders_decomposed_problem, +) +from .optimization import BlockBorderManagement, OptimizationProblem, build_problem from .output_values import OutputValues +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 new file mode 100644 index 00000000..e9dd9a6a --- /dev/null +++ b/src/andromede/simulation/benders_decomposed.py @@ -0,0 +1,219 @@ +# 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. + +""" +The xpansion module extends the optimization module +with Benders solver related functions +""" + +import json +import os +import pathlib +import subprocess +import sys +from typing import Any, Dict, List + +from andromede.simulation.optimization import ( + BlockBorderManagement, + OptimizationProblem, + build_problem, +) +from andromede.simulation.strategy import ( + InvestmentProblemStrategy, + OperationalProblemStrategy, +) +from andromede.simulation.time_block import TimeBlock +from andromede.study.data import DataBase +from andromede.study.network import Network +from andromede.utils import serialize + + +class BendersDecomposedProblem: + """ + A simpler interface for the Benders Decomposed problem + """ + + master: OptimizationProblem + subproblems: List[OptimizationProblem] + + def __init__( + self, master: OptimizationProblem, subproblems: List[OptimizationProblem] + ) -> None: + self.master = master + self.subproblems = subproblems + + def export_structure(self) -> str: + """ + Write the structure.txt file + """ + + if not self.subproblems: + # TODO For now, only one master and one subproblem + raise RuntimeError("Subproblem list must have at least one sub problem") + + # A mapping similar to the Xpansion mapping for keeping track of variable indexes + # in Master and Sub-problem files + problem_to_candidates: Dict[str, Dict[str, int]] = {} + candidates = set() + + problem_to_candidates["master"] = {} + for solver_var_info in self.master.context._solver_variables.values(): + if solver_var_info.is_in_objective: + problem_to_candidates["master"][ + solver_var_info.name + ] = solver_var_info.column_id + candidates.add(solver_var_info.name) + + for problem in self.subproblems: + problem_to_candidates[problem.name] = {} + + for solver_var_info in problem.context._solver_variables.values(): + if solver_var_info.name in candidates: + # If candidate was identified in master + problem_to_candidates[problem.name][ + solver_var_info.name + ] = solver_var_info.column_id + + structure_str = "" + for problem_name, candidate_to_index in problem_to_candidates.items(): + for candidate, index in candidate_to_index.items(): + structure_str += f"{problem_name:>50}{candidate:>50}{index:>10}\n" + + return structure_str + + def export_options( + self, *, solver_name: str = "XPRESS", log_level: int = 0 + ) -> Dict[str, Any]: + # Default values + options_values = { + "MAX_ITERATIONS": -1, + "ABSOLUTE_GAP": 1, + "RELATIVE_GAP": 1e-6, + "RELAXED_GAP": 1e-5, + "AGGREGATION": False, + "OUTPUTROOT": ".", + "TRACE": True, + "SLAVE_WEIGHT": "CONSTANT", + "SLAVE_WEIGHT_VALUE": 1, + "MASTER_NAME": "master", + "STRUCTURE_FILE": "structure.txt", + "INPUTROOT": ".", + "CSV_NAME": "benders_output_trace", + "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", + "MASTER_FORMULATION": "integer", + "SOLVER_NAME": solver_name, + "TIME_LIMIT": 1_000_000_000_000, + "LOG_LEVEL": log_level, + "LAST_MASTER_MPS": "master_last_iteration", + "LAST_MASTER_BASIS": "master_last_basis.bss", + } + return options_values + + 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( + "options.json", + json.dumps( + self.export_options(solver_name=solver_name, log_level=log_level), + indent=4, + ), + directory, + ) + + if is_debug: + serialize("master.lp", self.master.export_as_lp(), directory) + serialize("subproblem.lp", self.subproblems[0].export_as_lp(), directory) + + def run( + self, + *, + path: str = "outputs/lp", + solver_name: str = "XPRESS", + log_level: int = 0, + ) -> 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 + + 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) + + return res.returncode == 0 + + +def build_benders_decomposed_problem( + network: Network, + database: DataBase, + block: TimeBlock, + scenarios: int, + *, + border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, + solver_id: str = "GLOP", +) -> BendersDecomposedProblem: + """ + Entry point to build the xpansion problem for a time period + + Returns a Benders Decomposed problem + """ + + # Benders Decomposed Master Problem + master = build_problem( + network, + database, + block, + scenarios, + problem_name="master", + border_management=border_management, + solver_id=solver_id, + problem_strategy=InvestmentProblemStrategy(), + ) + + # 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(), + ) + + return BendersDecomposedProblem(master, [subproblem]) diff --git a/src/andromede/simulation/optimization.py b/src/andromede/simulation/optimization.py index 60d9916b..7e26dd38 100644 --- a/src/andromede/simulation/optimization.py +++ b/src/andromede/simulation/optimization.py @@ -18,7 +18,7 @@ from abc import ABC, abstractmethod from dataclasses import dataclass from enum import Enum -from typing import Dict, Iterable, List, Optional +from typing import Dict, Iterable, List, Optional, Type import ortools.linear_solver.pywraplp as lp @@ -40,6 +40,7 @@ from andromede.model.model import PortFieldId from andromede.simulation.linear_expression import LinearExpression, Term from andromede.simulation.linearize import linearize_expression +from andromede.simulation.strategy import MergedProblemStrategy, ModelSelectionStrategy from andromede.simulation.time_block import TimeBlock from andromede.study.data import DataBase from andromede.study.network import Component, Network @@ -222,10 +223,14 @@ def get_values(self, expression: ExpressionNode) -> TimestepValueProvider: ) def add_variable( - self, block_timestep: int, scenario: int, variable: lp.Variable + self, + block_timestep: int, + scenario: int, + model_var_name: str, + variable: lp.Variable, ) -> None: self.opt_context.register_component_variable( - block_timestep, scenario, self.component.id, variable.name(), variable + block_timestep, scenario, self.component.id, model_var_name, variable ) def get_variable( @@ -271,6 +276,20 @@ class BlockBorderManagement(Enum): CYCLE = "CYCLE" +@dataclass +class SolverVariableInfo: + """ + Helper class for constructing the structure file + for Benders solver. It keeps track of the corresponding + column of the variable in the MPS format as well as if it is + present in the objective function or not + """ + + name: str + column_id: int + is_in_objective: bool + + class OptimizationContext: """ Helper class to build the optimization problem. @@ -292,6 +311,7 @@ def __init__( self._scenarios = scenarios self._border_management = border_management self._component_variables: Dict[TimestepComponentVariableKey, lp.Variable] = {} + self._solver_variables: Dict[lp.Variable, SolverVariableInfo] = {} self._connection_fields_expressions: Dict[ PortFieldKey, List[ExpressionNode] ] = {} @@ -370,6 +390,10 @@ def register_component_variable( key = TimestepComponentVariableKey( component_id, variable_name, block_timestep, scenario ) + if key not in self._component_variables: + self._solver_variables[variable] = SolverVariableInfo( + variable.name(), len(self._solver_variables), False + ) self._component_variables[key] = variable def get_component_context(self, component: Component) -> ComponentContext: @@ -408,6 +432,23 @@ def _compute_indexing_structure( return constraint_indexing +def _instantiate_model_expression( + model_expression: ExpressionNode, + component_id: str, + optimization_context: OptimizationContext, +) -> ExpressionNode: + """ + Performs common operations that are necessary on model expressions before their actual use: + 1. add component ID for variables and parameters of THIS component + 2. replace port fields by their definition + """ + with_component = add_component_context(component_id, model_expression) + with_component_and_ports = resolve_port( + with_component, component_id, optimization_context.connection_fields_expressions + ) + return with_component_and_ports + + def _create_constraint( solver: lp.Solver, context: ComponentContext, @@ -449,6 +490,49 @@ def _create_constraint( ) +def _create_objective( + solver: lp.Solver, + opt_context: OptimizationContext, + component: Component, + component_context: ComponentContext, + objective_contribution: ExpressionNode, +) -> None: + instantiated_expr = _instantiate_model_expression( + objective_contribution, component.id, opt_context + ) + # We have already checked in the model creation that the objective contribution is neither indexed by time nor by scenario + linear_expr = component_context.linearize_expression(0, 0, instantiated_expr) + + obj: lp.Objective = solver.Objective() + for term in linear_expr.terms.values(): + # TODO : How to handle the scenario operator in a general manner ? + if isinstance(term.scenario_operator, Expectation): + weight = 1 / opt_context.scenarios + scenario_ids = range(opt_context.scenarios) + else: + weight = 1 + scenario_ids = range(1) + + for scenario in scenario_ids: + solver_vars = _get_solver_vars( + term, + opt_context, + 0, + scenario, + 0, + ) + + for solver_var in solver_vars: + opt_context._solver_variables[solver_var].is_in_objective = True + obj.SetCoefficient( + solver_var, + obj.GetCoefficient(solver_var) + weight * term.coefficient, + ) + + # This should have no effect on the optimization + obj.SetOffset(linear_expr.constant + obj.offset()) + + @dataclass class ConstraintData: name: str @@ -579,151 +663,147 @@ def make_constraint( return solver_constraints -def _create_variables( - network: Network, - opt_context: OptimizationContext, - solver: lp.Solver, -) -> None: - for component in network.all_components: - component_context = opt_context.get_component_context(component) - model = component.model +class OptimizationProblem: + name: str + solver: lp.Solver + context: OptimizationContext + strategy: ModelSelectionStrategy - for model_var in model.variables.values(): - var_indexing = IndexingStructure( - model_var.structure.time, model_var.structure.scenario - ) - instantiated_lb_expr = None - instantiated_ub_expr = None - if model_var.lower_bound: - instantiated_lb_expr = _instantiate_model_expression( - model_var.lower_bound, component.id, opt_context + def __init__( + self, + name: str, + solver: lp.Solver, + opt_context: OptimizationContext, + build_strategy: ModelSelectionStrategy = MergedProblemStrategy(), + ) -> None: + self.name = name + self.solver = solver + self.context = opt_context + self.strategy = build_strategy + + self._register_connection_fields_definitions() + self._create_variables() + self._create_constraints() + self._create_objectives() + + def _register_connection_fields_definitions(self) -> None: + for cnx in self.context.network.connections: + for field_name in list(cnx.master_port.keys()): + master_port = cnx.master_port[field_name] + port_definition = ( + master_port.component.model.port_fields_definitions.get( + PortFieldId( + port_name=master_port.port_id, field_name=field_name.name + ) + ) ) - if model_var.upper_bound: - instantiated_ub_expr = _instantiate_model_expression( - model_var.upper_bound, component.id, opt_context + expression_node = port_definition.definition # type: ignore + instantiated_expression = add_component_context( + master_port.component.id, expression_node + ) + self.context.register_connection_fields_expressions( + component_id=cnx.port1.component.id, + port_name=cnx.port1.port_id, + field_name=field_name.name, + expression=instantiated_expression, + ) + self.context.register_connection_fields_expressions( + component_id=cnx.port2.component.id, + port_name=cnx.port2.port_id, + field_name=field_name.name, + expression=instantiated_expression, ) - for block_timestep in opt_context.get_time_indices(var_indexing): - for scenario in opt_context.get_scenario_indices(var_indexing): - lower_bound = -solver.infinity() - upper_bound = solver.infinity() - if instantiated_lb_expr: - lower_bound = component_context.get_values( - instantiated_lb_expr - ).get_value(block_timestep, scenario) - if instantiated_ub_expr: - upper_bound = component_context.get_values( - instantiated_ub_expr - ).get_value(block_timestep, scenario) - - # TODO: Add BoolVar or IntVar if the variable is specified to be integer or bool - solver_var = solver.NumVar(lower_bound, upper_bound, model_var.name) - component_context.add_variable(block_timestep, scenario, solver_var) - - -def _register_connection_fields_definitions( - network: Network, opt_context: OptimizationContext -) -> None: - for cnx in network.connections: - for field_name in list(cnx.master_port.keys()): - master_port = cnx.master_port[field_name] - port_definition = master_port.component.model.port_fields_definitions.get( - PortFieldId(port_name=master_port.port_id, field_name=field_name.name) - ) - expression_node = port_definition.definition # type: ignore - instantiated_expression = add_component_context( - master_port.component.id, expression_node - ) - opt_context.register_connection_fields_expressions( - component_id=cnx.port1.component.id, - port_name=cnx.port1.port_id, - field_name=field_name.name, - expression=instantiated_expression, - ) - opt_context.register_connection_fields_expressions( - component_id=cnx.port2.component.id, - port_name=cnx.port2.port_id, - field_name=field_name.name, - expression=instantiated_expression, - ) + def _create_variables(self) -> None: + for component in self.context.network.all_components: + component_context = self.context.get_component_context(component) + model = component.model -def _create_constraints( - network: Network, opt_context: OptimizationContext, solver: lp.Solver -) -> None: - for component in network.all_components: - for constraint in component.model.get_all_constraints(): - instantiated_expr = _instantiate_model_expression( - constraint.expression, component.id, opt_context - ) - instantiated_lb = _instantiate_model_expression( - constraint.lower_bound, component.id, opt_context - ) - instantiated_ub = _instantiate_model_expression( - constraint.upper_bound, component.id, opt_context - ) - instantiated_constraint = Constraint( - name=constraint.name, - expression=instantiated_expr, - lower_bound=instantiated_lb, - upper_bound=instantiated_ub, - ) - _create_constraint( - solver, - opt_context.get_component_context(component), - instantiated_constraint, - ) + for model_var in self.strategy.get_variables(model): + var_indexing = IndexingStructure( + model_var.structure.time, model_var.structure.scenario + ) + instantiated_lb_expr = None + instantiated_ub_expr = None + if model_var.lower_bound: + instantiated_lb_expr = _instantiate_model_expression( + model_var.lower_bound, component.id, self.context + ) + if model_var.upper_bound: + instantiated_ub_expr = _instantiate_model_expression( + model_var.upper_bound, component.id, self.context + ) + for block_timestep in self.context.get_time_indices(var_indexing): + for scenario in self.context.get_scenario_indices(var_indexing): + lower_bound = -self.solver.infinity() + upper_bound = self.solver.infinity() + if instantiated_lb_expr: + lower_bound = component_context.get_values( + instantiated_lb_expr + ).get_value(block_timestep, scenario) + if instantiated_ub_expr: + upper_bound = component_context.get_values( + instantiated_ub_expr + ).get_value(block_timestep, scenario) + + # TODO: Add BoolVar or IntVar if the variable is specified to be integer or bool + # Externally, for the Solver, this variable will have a full name + # Internally, it will be indexed by a structure that into account + # the component id, variable name, timestep and scenario separately + solver_var = self.solver.NumVar( + lower_bound, + upper_bound, + f"{component.id}_{model_var.name}_{block_timestep}_{scenario}", + ) + component_context.add_variable( + block_timestep, scenario, model_var.name, solver_var + ) + def _create_constraints(self) -> None: + for component in self.context.network.all_components: + for constraint in self.strategy.get_constraints(component.model): + instantiated_expr = _instantiate_model_expression( + constraint.expression, component.id, self.context + ) + instantiated_lb = _instantiate_model_expression( + constraint.lower_bound, component.id, self.context + ) + instantiated_ub = _instantiate_model_expression( + constraint.upper_bound, component.id, self.context + ) -def _create_objective( - network: Network, - opt_context: OptimizationContext, - solver: lp.Solver, -) -> None: - for component in network.all_components: - component_context = opt_context.get_component_context(component) - model = component.model + instantiated_constraint = Constraint( + name=f"{component.id}_{constraint.name}", + expression=instantiated_expr, + lower_bound=instantiated_lb, + upper_bound=instantiated_ub, + ) + _create_constraint( + self.solver, + self.context.get_component_context(component), + instantiated_constraint, + ) - if model.objective_contribution is not None: - instantiated_expr = _instantiate_model_expression( - model.objective_contribution, component.id, opt_context - ) - # We have already checked in the model creation that the objective contribution is neither indexed by time nor by scenario - linear_expr = component_context.linearize_expression( - 0, 0, instantiated_expr - ) - obj: lp.Objective = solver.Objective() - for term in linear_expr.terms.values(): - # TODO : Comment gérer de manière générale les opérateurs sur les scénarios ? - if isinstance(term.scenario_operator, Expectation): - weight = 1 / opt_context.scenarios - scenario_ids = range(opt_context.scenarios) - else: - weight = 1 - scenario_ids = range(1) - - for scenario in scenario_ids: - solver_vars = _get_solver_vars( - term, - opt_context, - 0, - scenario, - 0, + def _create_objectives(self) -> None: + for component in self.context.network.all_components: + component_context = self.context.get_component_context(component) + model = component.model + + for objective in self.strategy.get_objectives(model): + if objective is not None: + _create_objective( + self.solver, + self.context, + component, + component_context, + objective, ) - for solver_var in solver_vars: - obj.SetCoefficient( - solver_var, - obj.GetCoefficient(solver_var) + weight * term.coefficient, - ) - - # This should have no effect on the optimization - obj.SetOffset(linear_expr.constant + obj.offset()) + def export_as_mps(self) -> str: + return self.solver.ExportModelAsMpsFormat(fixed_format=True, obfuscated=False) - -@dataclass(frozen=True) -class SolverAndContext: - solver: lp.Solver - context: "OptimizationContext" + def export_as_lp(self) -> str: + return self.solver.ExportModelAsLpFormat(obfuscated=False) def build_problem( @@ -731,9 +811,12 @@ def build_problem( database: DataBase, block: TimeBlock, scenarios: int, + *, + problem_name: str = "optimization_problem", border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, solver_id: str = "GLOP", -) -> SolverAndContext: + problem_strategy: ModelSelectionStrategy = MergedProblemStrategy(), +) -> OptimizationProblem: """ Entry point to build the optimization problem for a time period. """ @@ -745,25 +828,4 @@ def build_problem( network, database, block, scenarios, border_management ) - _register_connection_fields_definitions(network, opt_context) - _create_variables(network, opt_context, solver) - _create_constraints(network, opt_context, solver) - _create_objective(network, opt_context, solver) - return SolverAndContext(solver, opt_context) - - -def _instantiate_model_expression( - model_expression: ExpressionNode, - component_id: str, - optimization_context: OptimizationContext, -) -> ExpressionNode: - """ - Performs common operations that are necessary on model expressions before their actual use: - 1. add component ID for variables and parameters of THIS component - 2. replace port fields by their definition - """ - with_component = add_component_context(component_id, model_expression) - with_component_and_ports = resolve_port( - with_component, component_id, optimization_context.connection_fields_expressions - ) - return with_component_and_ports + return OptimizationProblem(problem_name, solver, opt_context, problem_strategy) diff --git a/src/andromede/simulation/output_values.py b/src/andromede/simulation/output_values.py index a569cdb2..d875bb8f 100644 --- a/src/andromede/simulation/output_values.py +++ b/src/andromede/simulation/output_values.py @@ -17,7 +17,7 @@ from dataclasses import dataclass, field from typing import Dict, List, Mapping, Optional, Tuple, TypeVar, Union, cast -from andromede.simulation.optimization import SolverAndContext +from andromede.simulation.optimization import OptimizationProblem from andromede.study.data import TimeScenarioIndex @@ -168,7 +168,7 @@ def var(self, variable_name: str) -> "OutputValues.Variable": self._variables[variable_name] = OutputValues.Variable(variable_name) return self._variables[variable_name] - problem: Optional[SolverAndContext] = field(default=None) + problem: Optional[OptimizationProblem] = field(default=None) _components: Dict[str, "OutputValues.Component"] = field( init=False, default_factory=dict ) diff --git a/src/andromede/simulation/strategy.py b/src/andromede/simulation/strategy.py new file mode 100644 index 00000000..75e34c65 --- /dev/null +++ b/src/andromede/simulation/strategy.py @@ -0,0 +1,82 @@ +# 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 typing import Generator, Optional + +from andromede.expression import ExpressionNode +from andromede.model import Constraint, Model, ProblemContext, Variable + + +class ModelSelectionStrategy(ABC): + """ + Abstract class to specify the strategy of the created problem. + Its derived classes select variables and constraints for the optimization problems: + - InvestmentProblemStrategy: Keep investment and coupling variables and constraints only for a BendersDecomposed master + - OperationalProblemStrategy: Keep operational and coupling variables and constraints only for a BendersDecomposed sub-problems + - MergedProblemStrategy: Keep all variables and constraints + """ + + def get_variables(self, model: Model) -> Generator[Variable, None, None]: + for variable in model.variables.values(): + if self._keep_from_context(variable.context): + yield variable + + def get_constraints(self, model: Model) -> Generator[Constraint, None, None]: + for constraint in model.get_all_constraints(): + if self._keep_from_context(constraint.context): + yield constraint + + @abstractmethod + def _keep_from_context(self, context: ProblemContext) -> bool: + ... + + @abstractmethod + def get_objectives( + self, model: Model + ) -> Generator[Optional[ExpressionNode], None, None]: + ... + + +class MergedProblemStrategy(ModelSelectionStrategy): + def _keep_from_context(self, context: ProblemContext) -> bool: + return True + + def get_objectives( + self, model: Model + ) -> Generator[Optional[ExpressionNode], None, None]: + yield model.objective_operational_contribution + yield model.objective_investment_contribution + + +class InvestmentProblemStrategy(ModelSelectionStrategy): + def _keep_from_context(self, context: ProblemContext) -> bool: + return ( + context == ProblemContext.INVESTMENT or context == ProblemContext.COUPLING + ) + + def get_objectives( + self, model: Model + ) -> Generator[Optional[ExpressionNode], None, None]: + yield model.objective_investment_contribution + + +class OperationalProblemStrategy(ModelSelectionStrategy): + def _keep_from_context(self, context: ProblemContext) -> bool: + return ( + context == ProblemContext.OPERATIONAL or context == ProblemContext.COUPLING + ) + + def get_objectives( + self, model: Model + ) -> Generator[Optional[ExpressionNode], None, None]: + yield model.objective_operational_contribution diff --git a/src/andromede/utils.py b/src/andromede/utils.py index 9a20a135..1b5c00b8 100644 --- a/src/andromede/utils.py +++ b/src/andromede/utils.py @@ -13,6 +13,7 @@ """ Module for technical utilities. """ +import pathlib from typing import Any, Callable, Dict, Generic, Optional, TypeVar T = TypeVar("T") @@ -41,3 +42,16 @@ def get_or_add(dictionary: Dict[K, V], key: K, default_factory: Supplier[V]) -> value = default_factory() dictionary[key] = value return value + + +def serialize(filename: str, message: str, path: pathlib.Path) -> None: + """ + Write message to path/filename + + Will throw an exception if it fails to create dir or ro open file + """ + path.mkdir(parents=True, exist_ok=True) + file = (path / filename).open(mode="w") + + with file: + file.write(message) diff --git a/tests/andromede/test_andromede.py b/tests/andromede/test_andromede.py index 5aeec26a..300ff88d 100644 --- a/tests/andromede/test_andromede.py +++ b/tests/andromede/test_andromede.py @@ -307,7 +307,9 @@ def test_variable_bound() -> None: definition=var("generation"), ) ], - objective_contribution=(param("cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("cost") * var("generation")) + .sum() + .expec(), ) network = create_one_node_network(generator_model) @@ -433,9 +435,13 @@ def test_min_up_down_times() -> None: PortRef(unsupplied_energy, "balance_port"), PortRef(node, "balance_port") ) - border_management = BlockBorderManagement.CYCLE - - problem = build_problem(network, database, time_block, scenarios, border_management) + problem = build_problem( + network, + database, + time_block, + scenarios, + border_management=BlockBorderManagement.CYCLE, + ) status = problem.solver.Solve() assert status == problem.solver.OPTIMAL @@ -511,9 +517,12 @@ def short_term_storage_base(efficiency: float, horizon: int) -> None: network.connect(PortRef(spillage, "balance_port"), PortRef(node, "balance_port")) network.connect(PortRef(unsupplied, "balance_port"), PortRef(node, "balance_port")) - border_management = BlockBorderManagement.CYCLE problem = build_problem( - network, database, time_blocks[0], scenarios, border_management + network, + database, + time_blocks[0], + scenarios, + border_management=BlockBorderManagement.CYCLE, ) status = problem.solver.Solve() diff --git a/tests/andromede/test_data.py b/tests/andromede/test_data.py index 5230c4c7..bd243c34 100644 --- a/tests/andromede/test_data.py +++ b/tests/andromede/test_data.py @@ -88,7 +88,9 @@ def mock_generator_with_fixed_scenario_time_varying_param() -> Model: name="Max generation", expression=var("generation") <= param("p_max") ) ], - objective_contribution=(param("cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("cost") * var("generation")) + .sum() + .expec(), ) return fixed_scenario_time_varying_param_generator @@ -114,7 +116,9 @@ def mock_generator_with_scenario_varying_fixed_time_param() -> Model: name="Max generation", expression=var("generation") <= param("p_max") ) ], - objective_contribution=(param("cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("cost") * var("generation")) + .sum() + .expec(), ) return scenario_varying_fixed_time_generator diff --git a/tests/andromede/test_electrolyzer.py b/tests/andromede/test_electrolyzer.py index 9a771da2..f6ea1d7c 100644 --- a/tests/andromede/test_electrolyzer.py +++ b/tests/andromede/test_electrolyzer.py @@ -63,7 +63,9 @@ definition=var("generation"), ) ], - objective_contribution=(param("cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("cost") * var("generation")) + .sum() + .expec(), ) H2_PORT = PortType(id="h2_port", fields=[PortField("flow")]) diff --git a/tests/andromede/test_model.py b/tests/andromede/test_model.py index 85238fbc..225e0451 100644 --- a/tests/andromede/test_model.py +++ b/tests/andromede/test_model.py @@ -188,7 +188,7 @@ def test_instantiating_a_model_with_non_linear_scenario_operator_in_the_objectiv _ = model( id="model_with_non_linear_op", variables=[float_variable("generation")], - objective_contribution=var("generation").variance(), + objective_operational_contribution=var("generation").variance(), ) assert str(exc.value) == "Objective contribution must be a linear expression." diff --git a/tests/andromede/test_output_values.py b/tests/andromede/test_output_values.py index feea3a1e..89b3637e 100644 --- a/tests/andromede/test_output_values.py +++ b/tests/andromede/test_output_values.py @@ -17,20 +17,18 @@ from andromede.simulation import OutputValues from andromede.simulation.optimization import ( OptimizationContext, - SolverAndContext, + OptimizationProblem, TimestepComponentVariableKey, ) def test_component_and_flow_output_object() -> None: mock_variable_component = Mock(spec=lp.Variable) - mock_variable_flow = Mock(spec=lp.Variable) + mock_problem = Mock(spec=OptimizationProblem) opt_context = Mock(spec=OptimizationContext) mock_variable_component.solution_value.side_effect = lambda: 1.0 - mock_variable_flow.solution_value.side_effect = lambda: -1.0 - opt_context.get_all_component_variables.return_value = { TimestepComponentVariableKey( component_id="component_id_test", @@ -48,8 +46,8 @@ def test_component_and_flow_output_object() -> None: opt_context.block_length.return_value = 1 - problem = SolverAndContext(mock_variable_flow, opt_context) - output = OutputValues(problem) + mock_problem.context = opt_context + output = OutputValues(mock_problem) test_output = OutputValues() assert output != test_output, f"Output is equal to empty output: {output}" diff --git a/tests/andromede/test_performance.py b/tests/andromede/test_performance.py index 25e6e2e1..f6ac6be1 100644 --- a/tests/andromede/test_performance.py +++ b/tests/andromede/test_performance.py @@ -52,7 +52,9 @@ def test_large_sum_with_loop() -> None: float_parameter(f"cost_{i}", IndexingStructure(False, False)) for i in range(1, nb_terms) ], - objective_contribution=sum([param(f"cost_{i}") for i in range(1, nb_terms)]), + objective_operational_contribution=sum( + [param(f"cost_{i}") for i in range(1, nb_terms)] + ), ) network = Network("test") @@ -86,7 +88,7 @@ def test_large_sum_outside_model() -> None: SIMPLE_COST_MODEL = model( id="SIMPLE_COST", parameters=[], - objective_contribution=literal(obj_coeff), + objective_operational_contribution=literal(obj_coeff), ) network = Network("test") diff --git a/tests/andromede/test_xpansion.py b/tests/andromede/test_xpansion.py index 5308a423..6e37faf3 100644 --- a/tests/andromede/test_xpansion.py +++ b/tests/andromede/test_xpansion.py @@ -20,17 +20,26 @@ DEMAND_MODEL, GENERATOR_MODEL, NODE_BALANCE_MODEL, + NODE_WITH_SPILL_AND_ENS_MODEL, ) from andromede.model import ( Constraint, Model, ModelPort, + ProblemContext, float_parameter, float_variable, + int_variable, model, ) from andromede.model.model import PortFieldDefinition, PortFieldId -from andromede.simulation import OutputValues, TimeBlock, build_problem +from andromede.simulation import ( + MergedProblemStrategy, + OutputValues, + TimeBlock, + build_benders_decomposed_problem, + build_problem, +) from andromede.study import ( Component, ConstantData, @@ -46,6 +55,10 @@ CONSTANT = IndexingStructure(False, False) FREE = IndexingStructure(True, True) +INVESTMENT = ProblemContext.INVESTMENT +OPERATIONAL = ProblemContext.OPERATIONAL +COUPLING = ProblemContext.COUPLING + @pytest.fixture def thermal_candidate() -> Model: @@ -62,6 +75,7 @@ def thermal_candidate() -> Model: lower_bound=literal(0), upper_bound=literal(1000), structure=CONSTANT, + context=COUPLING, ), ], ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")], @@ -76,12 +90,64 @@ def thermal_candidate() -> Model: name="Max generation", expression=var("generation") <= var("p_max") ) ], - objective_contribution=(param("invest_cost") * var("p_max")) - + (param("op_cost") * var("generation")).sum().expec(), + objective_operational_contribution=(param("op_cost") * var("generation")) + .sum() + .expec(), + objective_investment_contribution=param("invest_cost") * var("p_max"), ) return THERMAL_CANDIDATE +@pytest.fixture +def discrete_candidate() -> Model: + DISCRETE_CANDIDATE = model( + id="DISCRETE", + parameters=[ + float_parameter("op_cost", CONSTANT), + float_parameter("invest_cost", CONSTANT), + float_parameter("p_max_per_unit", CONSTANT), + ], + variables=[ + float_variable("generation", lower_bound=literal(0)), + float_variable( + "p_max", + lower_bound=literal(0), + structure=CONSTANT, + context=COUPLING, + ), + int_variable( + "nb_units", + lower_bound=literal(0), + upper_bound=literal(10), + structure=CONSTANT, + context=INVESTMENT, + ), + ], + ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")], + port_fields_definitions=[ + PortFieldDefinition( + port_field=PortFieldId("balance_port", "flow"), + definition=var("generation"), + ) + ], + constraints=[ + Constraint( + name="Max generation", expression=var("generation") <= var("p_max") + ), + Constraint( + name="Max investment", + expression=var("p_max") == param("p_max_per_unit") * var("nb_units"), + context=INVESTMENT, + ), + ], + objective_operational_contribution=(param("op_cost") * var("generation")) + .sum() + .expec(), + objective_investment_contribution=param("invest_cost") * var("p_max"), + ) + return DISCRETE_CANDIDATE + + @pytest.fixture def generator() -> Component: generator = create_component( @@ -97,6 +163,12 @@ def candidate(thermal_candidate: Model) -> Component: return candidate +@pytest.fixture +def cluster_candidate(discrete_candidate: Model) -> Component: + cluster = create_component(model=discrete_candidate, id="DISCRETE") + return cluster + + def test_generation_xpansion_single_time_step_single_scenario( generator: Component, candidate: Component, @@ -105,7 +177,7 @@ def test_generation_xpansion_single_time_step_single_scenario( Simple generation expansion problem on one node. One timestep, one scenario, one thermal cluster candidate. Demand = 300 - Generator : P_max : 100, Cost : 40 + Generator : P_max : 200, Cost : 40 Unsupplied energy : Cost : 1000 -> 100 of unsupplied energy @@ -144,33 +216,167 @@ def test_generation_xpansion_single_time_step_single_scenario( network.connect(PortRef(candidate, "balance_port"), PortRef(node, "balance_port")) scenarios = 1 + problem = build_problem( + network, + database, + TimeBlock(1, [0]), + scenarios, + problem_strategy=MergedProblemStrategy(), + ) + status = problem.solver.Solve() + + assert status == problem.solver.OPTIMAL + assert problem.solver.Objective().Value() == pytest.approx( + 490 * 100 + 100 * 10 + 200 * 40 + ) + + output = OutputValues(problem) + expected_output = OutputValues() + expected_output.component("G1").var("generation").value = 200.0 + expected_output.component("CAND").var("generation").value = 100.0 + expected_output.component("CAND").var("p_max").value = 100.0 + + assert output == expected_output, f"Output differs from expected: {output}" + + +def test_two_candidates_xpansion_single_time_step_single_scenario( + generator: Component, + candidate: Component, + cluster_candidate: Component, +) -> None: + """ + As before, 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. + + Demand = 400 + Generator : P_max : 200, Cost : 45 + Unsupplied energy : Cost : 1_000 + + -> 200 of unsupplied energy + -> Total cost without investment = 45 * 200 + 1_000 * 200 = 209_000 + + 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 (Single) + 200 * 100 (Cluster) -> investment + + 45 * 200 (Generator) + 10 * 100 (Single) + 10 * 100 (Cluster) -> operational + = 69_000 + 11_000 = 80_000 + """ + + database = DataBase() + database.add_data("D", "demand", ConstantData(400)) + + database.add_data("G1", "p_max", ConstantData(200)) + database.add_data("G1", "cost", ConstantData(45)) + + database.add_data("CAND", "op_cost", ConstantData(10)) + database.add_data("CAND", "invest_cost", ConstantData(490)) + + database.add_data("DISCRETE", "op_cost", ConstantData(10)) + database.add_data("DISCRETE", "invest_cost", ConstantData(200)) + database.add_data("DISCRETE", "p_max_per_unit", ConstantData(10)) + + demand = create_component(model=DEMAND_MODEL, id="D") + + node = Node(model=NODE_BALANCE_MODEL, id="N") + network = Network("test") + network.add_node(node) + network.add_component(demand) + network.add_component(generator) + network.add_component(candidate) + network.add_component(cluster_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")) + network.connect( + PortRef(cluster_candidate, "balance_port"), PortRef(node, "balance_port") + ) + scenarios = 1 + problem = build_problem(network, database, TimeBlock(1, [0]), scenarios) + status = problem.solver.Solve() assert status == problem.solver.OPTIMAL - assert problem.solver.Objective().Value() == 490 * 100 + 100 * 10 + 200 * 40 + assert problem.solver.Objective().Value() == pytest.approx( + (45 * 200) + (490 * 100 + 10 * 100) + (200 * 100 + 10 * 100) + ) output = OutputValues(problem) expected_output = OutputValues() expected_output.component("G1").var("generation").value = 200.0 expected_output.component("CAND").var("generation").value = 100.0 expected_output.component("CAND").var("p_max").value = 100.0 + expected_output.component("DISCRETE").var("generation").value = 100.0 + expected_output.component("DISCRETE").var("p_max").value = 100.0 + expected_output.component("DISCRETE").var("nb_units").value = 10.0 assert output == expected_output, f"Output differs from expected: {output}" +def test_model_export_xpansion_single_time_step_single_scenario( + 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 + """ + + database = DataBase() + database.add_data("D", "demand", ConstantData(400)) + + 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(45)) + + database.add_data("CAND", "op_cost", ConstantData(10)) + database.add_data("CAND", "invest_cost", ConstantData(490)) + + database.add_data("DISCRETE", "op_cost", ConstantData(10)) + database.add_data("DISCRETE", "invest_cost", ConstantData(200)) + database.add_data("DISCRETE", "p_max_per_unit", ConstantData(10)) + + 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.add_component(cluster_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")) + network.connect( + PortRef(cluster_candidate, "balance_port"), PortRef(node, "balance_port") + ) + scenarios = 1 + + xpansion = build_benders_decomposed_problem( + network, database, TimeBlock(1, [0]), scenarios + ) + assert xpansion.run() + + def test_generation_xpansion_two_time_steps_two_scenarios( generator: Component, candidate: Component, ) -> None: """ - Same as previous example but with two timesteps and two scenarios, in order to test the correct instanciation of the objective function + Same as previous example but with two timesteps and two scenarios, in order to test the correct instantiation of the objective function Demand = [300, 500] for scenario 1, [200, 400] for scenario 2 Generator : P_max : 200, Cost : 40 Unsupplied energy : Cost : 1000 - Scenarios équiprobables => Poids 0.5 + Equiprobable scenarios => Weights 0.5 -> 300 MW of unsupplied energy at step 2 for scenario 1 @@ -179,7 +385,7 @@ def test_generation_xpansion_two_time_steps_two_scenarios( Optimal investment : 300 MW -> Optimal cost = 490 * 300 (investment) - (poids 0.5 / scenario) + 0.5 * 10 * 300 (scenario 1, step 1) + (weights 0.5/scenario) + 0.5 * 10 * 300 (scenario 1, step 1) + 0.5 * (10 * 300 + 40 * 200) (scenario 1, step 2) + 0.5 * 10 * 200 (scenario 2, step 1) + 0.5 * (10 * 300 + 40 * 100) (scenario 2, step 2) @@ -227,9 +433,11 @@ def test_generation_xpansion_two_time_steps_two_scenarios( # assert ( # problem.solver.NumConstraints() == 3 * scenarios * horizon # ) # Flow balance, Max generation for each cluster - assert problem.solver.Objective().Value() == 490 * 300 + 0.5 * ( - 10 * 300 + 10 * 300 + 40 * 200 - ) + 0.5 * (10 * 200 + 10 * 300 + 40 * 100) + assert problem.solver.Objective().Value() == pytest.approx( + 490 * 300 + + 0.5 * (10 * 300 + 10 * 300 + 40 * 200) + + 0.5 * (10 * 200 + 10 * 300 + 40 * 100) + ) output = OutputValues(problem) expected_output = OutputValues()