From 5d760f5e3ec97f02e73dbf4dea5fb640f8dd730c Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Thu, 21 Mar 2024 14:08:54 +0100 Subject: [PATCH 01/14] Move pathway tests to functional folder --- tests/{andromede => functional}/test_generate_network_on_tree.py | 0 tests/{andromede => functional}/test_investment_pathway.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename tests/{andromede => functional}/test_generate_network_on_tree.py (100%) rename tests/{andromede => functional}/test_investment_pathway.py (100%) diff --git a/tests/andromede/test_generate_network_on_tree.py b/tests/functional/test_generate_network_on_tree.py similarity index 100% rename from tests/andromede/test_generate_network_on_tree.py rename to tests/functional/test_generate_network_on_tree.py diff --git a/tests/andromede/test_investment_pathway.py b/tests/functional/test_investment_pathway.py similarity index 100% rename from tests/andromede/test_investment_pathway.py rename to tests/functional/test_investment_pathway.py From fba71cd4a16b32f4bddd7597e8cf0b101e29f179 Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Thu, 21 Mar 2024 14:19:45 +0100 Subject: [PATCH 02/14] Add replicate method to classes to be able to copy them to tree --- src/andromede/expression/expression.py | 4 ++ src/andromede/model/constraint.py | 51 ++++++++++---------------- src/andromede/model/model.py | 20 +++++++--- src/andromede/model/parameter.py | 6 ++- src/andromede/model/variable.py | 7 +++- src/andromede/study/network.py | 48 ++++++++++++++++++------ tests/unittests/test_model.py | 24 ++++++------ 7 files changed, 96 insertions(+), 64 deletions(-) diff --git a/src/andromede/expression/expression.py b/src/andromede/expression/expression.py index ffbd4e22..4aa013f1 100644 --- a/src/andromede/expression/expression.py +++ b/src/andromede/expression/expression.py @@ -222,6 +222,10 @@ def literal(value: float) -> LiteralNode: return LiteralNode(value) +def is_unbound(expr: ExpressionNode) -> bool: + return isinstance(expr, LiteralNode) and (abs(expr.value) == float("inf")) + + @dataclass(frozen=True, eq=False) class UnaryOperatorNode(ExpressionNode): operand: ExpressionNode diff --git a/src/andromede/model/constraint.py b/src/andromede/model/constraint.py index 9e7c6050..ea87116c 100644 --- a/src/andromede/model/constraint.py +++ b/src/andromede/model/constraint.py @@ -10,19 +10,22 @@ # # This file is part of the Antares project. -from typing import Optional +from dataclasses import dataclass, field, replace +from typing import Any from andromede.expression.degree import is_constant from andromede.expression.expression import ( Comparator, ComparisonNode, ExpressionNode, + is_unbound, literal, ) from andromede.expression.print import print_expr from andromede.model.common import ProblemContext +@dataclass class Constraint: """ A constraint linking variables and parameters of a model together. @@ -32,55 +35,39 @@ class Constraint: name: str expression: ExpressionNode - lower_bound: ExpressionNode - upper_bound: ExpressionNode - context: ProblemContext + lower_bound: ExpressionNode = field(default=literal(-float("inf"))) + upper_bound: ExpressionNode = field(default=literal(float("inf"))) + context: ProblemContext = field(default=ProblemContext.OPERATIONAL) - def __init__( + def __post_init__( self, - name: str, - 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: + if isinstance(self.expression, ComparisonNode): + if not is_unbound(self.lower_bound) or not is_unbound(self.upper_bound): raise ValueError( "Both comparison between two expressions and a bound are specfied, set either only a comparison between expressions or a single linear expression with bounds." ) - merged_expr = expression.left - expression.right - self.expression = merged_expr - - if expression.comparator == Comparator.LESS_THAN: + if self.expression.comparator == Comparator.LESS_THAN: # lhs - rhs <= 0 self.upper_bound = literal(0) self.lower_bound = literal(-float("inf")) - elif expression.comparator == Comparator.GREATER_THAN: + elif self.expression.comparator == Comparator.GREATER_THAN: # lhs - rhs >= 0 self.lower_bound = literal(0) self.upper_bound = literal(float("inf")) else: # lhs - rhs == 0 self.lower_bound = literal(0) self.upper_bound = literal(0) + + self.expression = self.expression.left - self.expression.right + else: - for bound in [lower_bound, upper_bound]: - if bound is not None and not is_constant(bound): + for bound in [self.lower_bound, self.upper_bound]: + if not is_constant(bound): raise ValueError( f"The bounds of a constraint should not contain variables, {print_expr(bound)} was given." ) - self.expression = expression - if lower_bound is not None: - self.lower_bound = lower_bound - else: - self.lower_bound = literal(-float("inf")) - - if upper_bound is not None: - self.upper_bound = upper_bound - else: - self.upper_bound = literal(float("inf")) + def replicate(self, /, **changes: Any) -> "Constraint": + return replace(self, **changes) diff --git a/src/andromede/model/model.py b/src/andromede/model/model.py index 4a38e415..8871b2ab 100644 --- a/src/andromede/model/model.py +++ b/src/andromede/model/model.py @@ -16,11 +16,8 @@ defining parameters, variables, and equations. """ import itertools -from dataclasses import dataclass, field -from typing import Dict, Iterable, Optional - -from anytree import LevelOrderIter -from anytree import Node as TreeNode +from dataclasses import dataclass, field, replace +from typing import Any, Dict, Iterable, Optional from andromede.expression import ( AdditionNode, @@ -110,12 +107,18 @@ class ModelPort: port_type: PortType port_name: str + def replicate(self, /, **changes: Any) -> "ModelPort": + return replace(self, **changes) + @dataclass(frozen=True) class PortFieldId: port_name: str field_name: str + def replicate(self, /, **changes: Any) -> "PortFieldId": + return replace(self, **changes) + @dataclass(frozen=True) class PortFieldDefinition: @@ -129,6 +132,9 @@ class PortFieldDefinition: def __post_init__(self) -> None: _validate_port_field_expression(self) + def replicate(self, /, **changes: Any) -> "PortFieldDefinition": + return replace(self, **changes) + def port_field_def( port_name: str, field_name: str, definition: ExpressionNode @@ -186,6 +192,10 @@ def get_all_constraints(self) -> Iterable[Constraint]: self.binding_constraints.values(), self.constraints.values() ) + def replicate(self, /, **changes: Any) -> "Model": + # Shallow copy + return replace(self, **changes) + def model( id: str, diff --git a/src/andromede/model/parameter.py b/src/andromede/model/parameter.py index 745d5a40..fdba78c4 100644 --- a/src/andromede/model/parameter.py +++ b/src/andromede/model/parameter.py @@ -10,7 +10,8 @@ # # This file is part of the Antares project. -from dataclasses import dataclass +from dataclasses import dataclass, replace +from typing import Any from andromede.expression.indexing_structure import IndexingStructure from andromede.model.common import ValueType @@ -28,6 +29,9 @@ class Parameter: type: ValueType structure: IndexingStructure + def replicate(self, /, **changes: Any) -> "Parameter": + return replace(self, **changes) + def int_parameter( name: str, diff --git a/src/andromede/model/variable.py b/src/andromede/model/variable.py index 7f65b2c4..ca4ca0cd 100644 --- a/src/andromede/model/variable.py +++ b/src/andromede/model/variable.py @@ -10,8 +10,8 @@ # # This file is part of the Antares project. -from dataclasses import dataclass -from typing import Optional +from dataclasses import dataclass, replace +from typing import Any, Optional from andromede.expression import ExpressionNode from andromede.expression.degree import is_constant @@ -38,6 +38,9 @@ def __post_init__(self) -> None: if self.upper_bound and not is_constant(self.upper_bound): raise ValueError("Upper bounds of variables must be constant") + def replicate(self, /, **changes: Any) -> "Variable": + return replace(self, **changes) + def int_variable( name: str, diff --git a/src/andromede/study/network.py b/src/andromede/study/network.py index 44b77974..4bf0c5d1 100644 --- a/src/andromede/study/network.py +++ b/src/andromede/study/network.py @@ -15,8 +15,8 @@ including nodes, links, and components (model instantations). """ import itertools -from dataclasses import dataclass -from typing import Dict, Iterable, List +from dataclasses import dataclass, field, replace +from typing import Any, Dict, Iterable, List, cast from andromede.model import PortField, PortType from andromede.model.model import Model, PortFieldId @@ -32,6 +32,9 @@ class Component: model: Model id: str + def replicate(self, /, **changes: Any) -> "Component": + return replace(self, **changes) + def create_component(model: Model, id: str) -> Component: return Component(model=model, id=id) @@ -46,6 +49,10 @@ class Node(Component): pass +def create_node(model: Model, id: str) -> Node: + return Node(model=model, id=id) + + @dataclass(frozen=True) class PortRef: component: Component @@ -56,12 +63,9 @@ class PortRef: class PortsConnection: port1: PortRef port2: PortRef - master_port: Dict[PortField, PortRef] + master_port: Dict[PortField, PortRef] = field(init=False, default_factory=dict) - def __init__(self, port1: PortRef, port2: PortRef): - self.port1 = port1 - self.port2 = port2 - self.master_port = {} + def __post_init__(self) -> None: self.__validate_ports() def __validate_ports(self) -> None: @@ -106,6 +110,10 @@ def get_port_type(self) -> PortType: raise ValueError(f"Missing port: {port_1}") return port_1.port_type + def replicate(self, /, **changes: Any) -> "PortsConnection": + # Shallow copy + return replace(self, **changes) + @dataclass class Network: @@ -113,11 +121,10 @@ class Network: Network model: simply nodes, links, and components. """ - def __init__(self, id: str): - self.id: str = id - self._nodes: Dict[str, Node] = {} - self._components: Dict[str, Component] = {} - self._connections: List[PortsConnection] = [] + id: str + _nodes: Dict[str, Node] = field(init=False, default_factory=dict) + _components: Dict[str, Component] = field(init=False, default_factory=dict) + _connections: List[PortsConnection] = field(init=False, default_factory=list) def _check_node_exists(self, node_id: str) -> None: if node_id not in self._nodes: @@ -162,3 +169,20 @@ def connect(self, port1: PortRef, port2: PortRef) -> None: @property def connections(self) -> Iterable[PortsConnection]: return self._connections + + def get_connection(self, idx: int) -> PortsConnection: + return self._connections[idx] + + def replicate(self, /, **changes: Any) -> "Network": + replica = replace(self, **changes) + + for node in self.nodes: + replica.add_node(cast(Node, node.replicate())) + + for component in self.components: + replica.add_component(component.replicate()) + + for connection in self.connections: + replica._connections.append(connection.replicate()) + + return replica diff --git a/tests/unittests/test_model.py b/tests/unittests/test_model.py index 225e0451..b5efe767 100644 --- a/tests/unittests/test_model.py +++ b/tests/unittests/test_model.py @@ -42,8 +42,8 @@ ( "my_constraint", 2 * var("my_var"), - None, - None, + literal(-float("inf")), + literal(float("inf")), "my_constraint", 2 * var("my_var"), literal(-float("inf")), @@ -52,8 +52,8 @@ ( "my_constraint", 2 * var("my_var") <= param("p"), - None, - None, + literal(-float("inf")), + literal(float("inf")), "my_constraint", 2 * var("my_var") - param("p"), literal(-float("inf")), @@ -62,8 +62,8 @@ ( "my_constraint", 2 * var("my_var") >= param("p"), - None, - None, + literal(-float("inf")), + literal(float("inf")), "my_constraint", 2 * var("my_var") - param("p"), literal(0), @@ -72,8 +72,8 @@ ( "my_constraint", 2 * var("my_var") == param("p"), - None, - None, + literal(-float("inf")), + literal(float("inf")), "my_constraint", 2 * var("my_var") - param("p"), literal(0), @@ -82,8 +82,8 @@ ( "my_constraint", 2 * var("my_var").expec() == param("p"), - None, - None, + literal(-float("inf")), + literal(float("inf")), "my_constraint", 2 * var("my_var").expec() - param("p"), literal(0), @@ -92,8 +92,8 @@ ( "my_constraint", 2 * var("my_var").shift(-1) == param("p"), - None, - None, + literal(-float("inf")), + literal(float("inf")), "my_constraint", 2 * var("my_var").shift(-1) - param("p"), literal(0), From 3e92e2025f743c0563401b2789af372c4c4b7385 Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Thu, 21 Mar 2024 17:27:19 +0100 Subject: [PATCH 03/14] Remove integration test from functional xpansion test --- tests/functional/test_xpansion.py | 60 ------------------------------- 1 file changed, 60 deletions(-) diff --git a/tests/functional/test_xpansion.py b/tests/functional/test_xpansion.py index d02b4e74..364a1a05 100644 --- a/tests/functional/test_xpansion.py +++ b/tests/functional/test_xpansion.py @@ -20,7 +20,6 @@ DEMAND_MODEL, GENERATOR_MODEL, NODE_BALANCE_MODEL, - NODE_WITH_SPILL_AND_ENS, THERMAL_CANDIDATE, ) from andromede.model import ( @@ -40,11 +39,6 @@ TimeBlock, build_problem, ) -from andromede.simulation.benders_decomposed import build_benders_decomposed_problem -from andromede.simulation.decision_tree import ( - create_network_on_tree, - create_single_node_decision_tree, -) from andromede.study import ( Component, ConstantData, @@ -286,60 +280,6 @@ def test_two_candidates_xpansion_single_time_step_single_scenario( 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, - demand: 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("CAND", "max_invest", ConstantData(1000)) - - 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)) - - node = Node(model=NODE_WITH_SPILL_AND_ENS, 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") - ) - - blocks = [TimeBlock(1, [0])] - scenarios = 1 - - configured_tree = create_single_node_decision_tree(blocks, scenarios) - tree_node_to_network = create_network_on_tree(network, configured_tree.root) - - xpansion = build_benders_decomposed_problem( - tree_node_to_network, database, configured_tree - ) - assert xpansion.run() - - def test_generation_xpansion_two_time_steps_two_scenarios( generator: Component, candidate: Component, From d12d77f3e29fac08b5c0b8fdb68b2685d10a1ddb Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Thu, 21 Mar 2024 17:31:06 +0100 Subject: [PATCH 04/14] Modify DecisionTree --- .../simulation/benders_decomposed.py | 20 +- src/andromede/simulation/decision_tree.py | 176 ++++-------------- .../test_generate_network_on_tree.py | 39 ++-- tests/functional/test_investment_pathway.py | 36 ++-- tests/integration/test_benders_decomposed.py | 20 +- 5 files changed, 89 insertions(+), 202 deletions(-) diff --git a/src/andromede/simulation/benders_decomposed.py b/src/andromede/simulation/benders_decomposed.py index e94d968e..4429a1b5 100644 --- a/src/andromede/simulation/benders_decomposed.py +++ b/src/andromede/simulation/benders_decomposed.py @@ -21,7 +21,7 @@ from anytree import Node as TreeNode from andromede.model.model import Model -from andromede.simulation.decision_tree import ConfiguredTree, create_master_network +from andromede.simulation.decision_tree import DecisionTreeNode, create_master_network from andromede.simulation.optimization import ( BlockBorderManagement, OptimizationProblem, @@ -208,9 +208,8 @@ def run( def build_benders_decomposed_problem( - network_on_tree: Dict[TreeNode, Network], + decision_tree_root: DecisionTreeNode, database: DataBase, - configured_tree: ConfiguredTree, *, decision_coupling_model: Optional[Model] = None, border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, @@ -222,7 +221,7 @@ def build_benders_decomposed_problem( Returns a Benders Decomposed problem """ - master_network = create_master_network(network_on_tree, decision_coupling_model) + master_network = create_master_network(decision_tree_root, decision_coupling_model) # Benders Decomposed Master Problem master = build_problem( @@ -240,19 +239,16 @@ def build_benders_decomposed_problem( # Benders Decomposed Sub-problems subproblems = [] - for ( - tree_node, - time_scenario_config, - ) in configured_tree.node_to_config.items(): - for block in time_scenario_config.blocks: + for tree_node in decision_tree_root.traverse(): + for block in tree_node.config.blocks: # Xpansion Sub-problems subproblems.append( build_problem( - network_on_tree[tree_node], + tree_node.network, database, block, - time_scenario_config.scenarios, - problem_name=f"subproblem_{tree_node.name}_{block.id}", + tree_node.config.scenarios, + problem_name=f"subproblem_{tree_node.id}_{block.id}", solver_id=solver_id, problem_strategy=OperationalProblemStrategy(), ) diff --git a/src/andromede/simulation/decision_tree.py b/src/andromede/simulation/decision_tree.py index 267c71ee..8c883910 100644 --- a/src/andromede/simulation/decision_tree.py +++ b/src/andromede/simulation/decision_tree.py @@ -10,19 +10,14 @@ # # This file is part of the Antares project. -import dataclasses -from dataclasses import dataclass, field -from typing import Dict, Iterable, List, Optional +from dataclasses import dataclass +from typing import Generator, Iterable, List, Optional -from anytree import LevelOrderIter -from anytree import Node as TreeNode +from anytree import LevelOrderIter, NodeMixin -from andromede.expression.expression import ExpressionNode -from andromede.model.constraint import Constraint -from andromede.model.model import Model, PortFieldDefinition, PortFieldId, model -from andromede.model.variable import Variable +from andromede.model.model import Model from andromede.simulation.time_block import TimeBlock -from andromede.study.network import Component, Network, Node, create_component +from andromede.study.network import Network @dataclass(frozen=True) @@ -31,142 +26,47 @@ class InterDecisionTimeScenarioConfig: scenarios: int -@dataclass(frozen=True) -class ConfiguredTree: - node_to_config: Dict[TreeNode, InterDecisionTimeScenarioConfig] - root: TreeNode = field(init=False) - - def __post_init__(self) -> None: - # Stores the root, by getting it from any tree node - object.__setattr__(self, "root", next(iter(self.node_to_config.keys())).root) - - -def create_single_node_decision_tree( - blocks: List[TimeBlock], scenarios: int -) -> ConfiguredTree: - time_scenario_config = InterDecisionTimeScenarioConfig(blocks, scenarios) - - root = TreeNode("root") - configured_tree = ConfiguredTree( - { - root: time_scenario_config, - }, - ) - - return configured_tree - - -def _generate_tree_variables( - variables: Dict[str, Variable], tree_node: TreeNode -) -> Iterable[Variable]: - tree_variables = [] - for variable in variables.values(): - # Works as we do not allow variables in bounds, hence no problem to copy the corresponding expression nodes as is. If we had variables, we would have to replace the variable names by the ones with tree node information. - tree_variables.append( - dataclasses.replace(variable, name=f"{tree_node.name}_{variable.name}") - ) - return tree_variables - - -def _generate_tree_constraints( - constraints: Dict[str, Constraint], tree_node: TreeNode -) -> Iterable[Constraint]: - # Goal is to replace variables in constraint, lower bound and upper bound with node variable - raise NotImplementedError() - +class DecisionTreeNode(NodeMixin): + id: str + config: InterDecisionTimeScenarioConfig + network: Network -def _generate_tree_expression( - expression: Optional[ExpressionNode], tree_node: TreeNode -) -> ExpressionNode: - # Goal is to replace variables with node variable - # Create a copy visitor to do so - raise NotImplementedError() + def __init__( + self, + id: str, + config: InterDecisionTimeScenarioConfig, + network: Network = Network(""), + parent: Optional["DecisionTreeNode"] = None, + children: Optional[Iterable["DecisionTreeNode"]] = None, + ) -> None: + self.id = id + self.config = config + self.network = network + self.parent = parent + if children: + self.children = children + def traverse(self) -> Generator["DecisionTreeNode", None, None]: + yield from LevelOrderIter(self) -def _generate_tree_port_field_definition( - port_field_definition: Dict[PortFieldId, PortFieldDefinition], tree_node: TreeNode -) -> Iterable[PortFieldDefinition]: - # Goal is to replace variables in the expression defining the port by node variable - raise NotImplementedError() +def replicate_network_from_root(root: DecisionTreeNode) -> None: + if root.size == 1: + # Nothing to replicate. Just past down the network + return -def _generate_tree_model( - tree_node: TreeNode, - component: Component, -) -> Model: - variables = _generate_tree_variables( - component.model.variables, - tree_node, - ) - constraints = _generate_tree_constraints(component.model.constraints, tree_node) - binding_constraints = _generate_tree_constraints( - component.model.binding_constraints, tree_node - ) - objective_operational_contribution = _generate_tree_expression( - component.model.objective_operational_contribution, tree_node - ) - objective_investment_contribution = _generate_tree_expression( - component.model.objective_investment_contribution, tree_node - ) - port_fields_definitions = _generate_tree_port_field_definition( - component.model.port_fields_definitions, tree_node - ) - tree_model = model( - id=f"{tree_node.name}_{component.model.id}", - constraints=constraints, - binding_constraints=binding_constraints, - parameters=component.model.parameters.values(), - variables=variables, - objective_operational_contribution=objective_operational_contribution, - objective_investment_contribution=objective_investment_contribution, - inter_block_dyn=component.model.inter_block_dyn, - ports=component.model.ports.values(), - port_fields_definitions=port_fields_definitions, - ) - - return tree_model - - -def _generate_network_on_node(network: Network, tree_node: TreeNode) -> Network: - tree_node_network = Network(tree_node.name) - - for component in network.all_components: - tree_node_model = _generate_tree_model(tree_node, component) - - # It would be nice to have the same treatment for nodes and components as they are actually the same thing... - if isinstance(component, Node): - network_node = Node(tree_node_model, id=f"{tree_node.name}_{component.id}") - tree_node_network.add_node(network_node) - else: - tree_node_component = create_component( - tree_node_model, id=f"{tree_node.name}_{component.id}" - ) - tree_node_network.add_component(tree_node_component) - - for connection in network.connections: - tree_node_network.connect(connection.port1, connection.port2) - return tree_node_network - - -def create_network_on_tree(network: Network, tree: TreeNode) -> Dict[TreeNode, Network]: - # On crée un gros modèle en dupliquant les variables; contraintes, etc à chaque noeud de l'arbre. - # Pour le master on peut : - # - Utiliser uniquement les variables, contraintes, etc dont on va avoir besoin dans la construction du problème -> nécessite déjà d'avoir des infos sur la construction des problèmes alors qu'on agit au niveau modèle ici - # - Dupliquer tout le modèle, permet de mutualiser du code avec la partie composant par noeud et plus lisible. Seul inconvénient, modèle master un peu trop riche, pas besoin des infos "opérationnelles". Mais les modèles ne sont pas très "lourds" donc on peut se le permettre. C'est l'option choisie ici. - if tree.size == 1: - return {tree: network} else: - node_to_network = {} - for tree_node in LevelOrderIter(tree): - node_to_network[tree_node] = _generate_network_on_node(network, tree_node) - return node_to_network + # Replicates the network and updates their id to take into account the tree node's id + original_network_id = root.network.id + for tree_node in LevelOrderIter(root): + tree_node.network = root.network.replicate( + id=f"{tree_node.id}_{original_network_id}" + ) def create_master_network( - tree_node_to_network: Dict[TreeNode, Network], + root: DecisionTreeNode, decision_coupling_model: Optional[Model], ) -> Network: - # Current implementation so that tests pass for trees with one investment nodes (in test_xpansion) - # The final implementation should gather all networks from tree nodes and connect the models with the decision coupling model (with ports) - root = next(iter(tree_node_to_network.keys())).root - return tree_node_to_network[root] + # TODO + return root.network diff --git a/tests/functional/test_generate_network_on_tree.py b/tests/functional/test_generate_network_on_tree.py index b2e516a6..c77111d3 100644 --- a/tests/functional/test_generate_network_on_tree.py +++ b/tests/functional/test_generate_network_on_tree.py @@ -11,25 +11,34 @@ # This file is part of the Antares project. -from anytree import Node as TreeNode - -from andromede.libs.standard import THERMAL_CLUSTER_MODEL_HD -from andromede.simulation.decision_tree import _generate_tree_model -from andromede.study.network import create_component +from andromede.simulation import TimeBlock +from andromede.simulation.decision_tree import ( + DecisionTreeNode, + InterDecisionTimeScenarioConfig, + replicate_network_from_root, +) +from andromede.study.network import Network def test_generate_model_on_node() -> None: - thermal = create_component(model=THERMAL_CLUSTER_MODEL_HD, id="thermal") + scenarios = 1 + blocks = [TimeBlock(1, [0])] + config = InterDecisionTimeScenarioConfig(blocks, scenarios) + + network = Network("network_id") + tree_root = DecisionTreeNode("root", config, network) + + assert tree_root.id == "root" + assert tree_root.parent is None + assert not tree_root.children # No children + + child = DecisionTreeNode("child", config, parent=tree_root) - tree_node_id = "2030" - tree_node_model = _generate_tree_model(TreeNode(tree_node_id), thermal) + print(child.parent) - # How to compare model efficiently with only change in name ? - assert tree_node_model.id == f"{tree_node_id}_{thermal.id}" + assert child.parent == tree_root + assert child in tree_root.children - for variable in thermal.model.variables.values(): - assert f"{tree_node_id}_{variable.name}" in tree_node_model.variables + replicate_network_from_root(tree_root) - # Create dedicated function - tree_variable = tree_node_model.variables[f"{tree_node_id}_{variable.name}"] - # assert + assert child.network.id == "child_network_id" diff --git a/tests/functional/test_investment_pathway.py b/tests/functional/test_investment_pathway.py index 1b7b0052..33d4f703 100644 --- a/tests/functional/test_investment_pathway.py +++ b/tests/functional/test_investment_pathway.py @@ -22,9 +22,9 @@ from andromede.model.model import model from andromede.simulation.benders_decomposed import build_benders_decomposed_problem from andromede.simulation.decision_tree import ( - ConfiguredTree, + DecisionTreeNode, InterDecisionTimeScenarioConfig, - create_network_on_tree, + replicate_network_from_root, ) from andromede.simulation.time_block import TimeBlock from andromede.study.data import ConstantData, DataBase, TreeData @@ -176,11 +176,6 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( ) database.add_data("BASE", "cost", ConstantData(5)) - # Fonction qui crée les composants / noeud en fonction de l'arbre et du Database initial / modèles + générer les contraintes couplantes temporelles trajectoire + actualisation + - # contraintes industrielles liées à l'arbre ? - # Test mode peigne - # Générer le modèle "couplant" - network = Network("test") network.add_node(node) network.add_component(demand) @@ -195,26 +190,17 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( [TimeBlock(0, [0])], scenarios ) - root = TreeNode("2030") - new_base = TreeNode("2040_new_base", parent=root) - no_base = TreeNode("2040_no_base", parent=root) - configured_tree = ConfiguredTree( - { - root: time_scenario_config, - new_base: time_scenario_config, - no_base: time_scenario_config, - }, + decision_tree_root = DecisionTreeNode("2030", time_scenario_config, network) + decision_tree_new_base = DecisionTreeNode( + "2040_new_base", time_scenario_config, parent=decision_tree_root + ) + decision_tree_no_base = DecisionTreeNode( + "2040_no_base", time_scenario_config, parent=decision_tree_root ) + replicate_network_from_root(decision_tree_root) decision_coupling_model = model("DECISION_COUPLING") - tree_node_to_network = create_network_on_tree(network, configured_tree.root) - - problems = build_benders_decomposed_problem( - tree_node_to_network, - database, - configured_tree, - decision_coupling_model=decision_coupling_model, + xpansion = build_benders_decomposed_problem( + decision_tree_root, database, decision_coupling_model=decision_coupling_model ) - - # Réfléchir à la représentation des variables dans l'arbre diff --git a/tests/integration/test_benders_decomposed.py b/tests/integration/test_benders_decomposed.py index a6011161..a2ef3269 100644 --- a/tests/integration/test_benders_decomposed.py +++ b/tests/integration/test_benders_decomposed.py @@ -38,8 +38,8 @@ build_benders_decomposed_problem, ) from andromede.simulation.decision_tree import ( - create_network_on_tree, - create_single_node_decision_tree, + DecisionTreeNode, + InterDecisionTimeScenarioConfig, ) from andromede.study import ( Component, @@ -233,12 +233,10 @@ def test_benders_decomposed_integration( scenarios = 1 blocks = [TimeBlock(1, [0])] - configured_tree = create_single_node_decision_tree(blocks, scenarios) - tree_node_to_network = create_network_on_tree(network, configured_tree.root) + config = InterDecisionTimeScenarioConfig(blocks, scenarios) + decision_tree_root = DecisionTreeNode("root", config, network) - xpansion = build_benders_decomposed_problem( - tree_node_to_network, database, configured_tree - ) + xpansion = build_benders_decomposed_problem(decision_tree_root, database) data = { "solution": { @@ -325,12 +323,10 @@ def test_benders_decomposed_multi_time_block_single_scenario( scenarios = 1 blocks = [TimeBlock(1, [0]), TimeBlock(2, [1])] - configured_tree = create_single_node_decision_tree(blocks, scenarios) - tree_node_to_network = create_network_on_tree(network, configured_tree.root) + config = InterDecisionTimeScenarioConfig(blocks, scenarios) + decision_tree_root = DecisionTreeNode("root", config, network) - xpansion = build_benders_decomposed_problem( - tree_node_to_network, database, configured_tree - ) + xpansion = build_benders_decomposed_problem(decision_tree_root, database) data = { "solution": { From c42c48ff63f213ba74e79371f645bd57c23ee60f Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Thu, 21 Mar 2024 17:32:26 +0100 Subject: [PATCH 05/14] Move decision_tree test to unit tests folder --- tests/{functional => unittests}/test_generate_network_on_tree.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tests/{functional => unittests}/test_generate_network_on_tree.py (100%) diff --git a/tests/functional/test_generate_network_on_tree.py b/tests/unittests/test_generate_network_on_tree.py similarity index 100% rename from tests/functional/test_generate_network_on_tree.py rename to tests/unittests/test_generate_network_on_tree.py From 35b67584aa7b76bd3688bc38d2ad7419804e5b83 Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Fri, 22 Mar 2024 17:26:32 +0100 Subject: [PATCH 06/14] Ideas for the coupling_model --- src/andromede/simulation/decision_tree.py | 50 ++++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/src/andromede/simulation/decision_tree.py b/src/andromede/simulation/decision_tree.py index 8c883910..5737ade5 100644 --- a/src/andromede/simulation/decision_tree.py +++ b/src/andromede/simulation/decision_tree.py @@ -16,6 +16,7 @@ from anytree import LevelOrderIter, NodeMixin from andromede.model.model import Model +from andromede.simulation.strategy import InvestmentProblemStrategy from andromede.simulation.time_block import TimeBlock from andromede.study.network import Network @@ -68,5 +69,52 @@ def create_master_network( root: DecisionTreeNode, decision_coupling_model: Optional[Model], ) -> Network: - # TODO + # TODO Use ports for coupling different models across the decision tree + """ + Each candidate model should have one of these ports. + As for balance, they all should have a common name like "investment" and they could be defined as: + + ports=[ModelPort(port_type=INVESTMENT_PORT_TYPE, port_name="investment_port")], + port_fields_definitions=[ + PortFieldDefinition( + port_field=PortFieldId("investment_port", "investment"), + definition=param("p_max"), ---> p_max if we are investing the p_max for instance + ) + ], + + Then, here, we would traverse the tree and create a Network that only keeps + variables allowed in the InvestmentProblemStrategy (we would need to change their names to something like + p_max_root, p_max_child_id, p_max_child_id2, ...) + We would connect them by some constraints using these ports to the decision_coupling_model (a parameter here), + so we could have something like: + + master_network.connect(PortRef(candidate_on_root, "investment_port"), PortRef(coupling_model, "investment_port")) + master_network.connect(PortRef(candidate_on_child, "investment_port"), PortRef(coupling_model, "investment_port")) + + On the coupling model, we would have something like for the nodes: + ports=[ModelPort(port_type=INVESTMENT_PORT_TYPE, port_name="investment_port")], + binding_constraints=[ + Constraint( + name="Pathway_Investment", + expression=port_field("investment_port", "flow").some_operator(), + ) + ], + + Maybe we would have to define a operator to represent consecutive inequalities ? + Or create one coupling model per pair parent-child so the expression would become + port_field("investment_port", "flow").sum_connections() <= 0 ? In this case, we would have to define + a negative and a positive value + + To resume, a network here would need: + - All candidates on all tree nodes (so we don't make great changes on build_problem) + - Update investment variable names to show their corresponding tree node; + - Add a investment port to each model + - A coupling model with the good ports to connect to + - The binding constraint that ties everything + """ + + strategy = InvestmentProblemStrategy() + for tree_node in root.traverse(): + ... + return root.network From 393d39a2333bcd92228557452890bedbb2701ee2 Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Wed, 17 Apr 2024 14:01:21 +0200 Subject: [PATCH 07/14] Write subproblem with tree --- .../simulation/benders_decomposed.py | 1 + src/andromede/simulation/decision_tree.py | 21 +++++--- src/andromede/simulation/optimization.py | 26 +++------- src/andromede/study/data.py | 50 ++++++------------- tests/functional/test_investment_pathway.py | 42 +++++++++++++--- tests/integration/test_benders_decomposed.py | 2 +- .../test_generate_network_on_tree.py | 24 +++++---- 7 files changed, 87 insertions(+), 79 deletions(-) diff --git a/src/andromede/simulation/benders_decomposed.py b/src/andromede/simulation/benders_decomposed.py index 4429a1b5..5e42306f 100644 --- a/src/andromede/simulation/benders_decomposed.py +++ b/src/andromede/simulation/benders_decomposed.py @@ -251,6 +251,7 @@ def build_benders_decomposed_problem( problem_name=f"subproblem_{tree_node.id}_{block.id}", solver_id=solver_id, problem_strategy=OperationalProblemStrategy(), + decision_tree_node=tree_node.id, ) ) diff --git a/src/andromede/simulation/decision_tree.py b/src/andromede/simulation/decision_tree.py index 5737ade5..505a6904 100644 --- a/src/andromede/simulation/decision_tree.py +++ b/src/andromede/simulation/decision_tree.py @@ -51,7 +51,7 @@ def traverse(self) -> Generator["DecisionTreeNode", None, None]: yield from LevelOrderIter(self) -def replicate_network_from_root(root: DecisionTreeNode) -> None: +def replicate_network_from(root: DecisionTreeNode) -> None: if root.size == 1: # Nothing to replicate. Just past down the network return @@ -59,7 +59,7 @@ def replicate_network_from_root(root: DecisionTreeNode) -> None: else: # Replicates the network and updates their id to take into account the tree node's id original_network_id = root.network.id - for tree_node in LevelOrderIter(root): + for tree_node in root.traverse(): tree_node.network = root.network.replicate( id=f"{tree_node.id}_{original_network_id}" ) @@ -112,9 +112,18 @@ def create_master_network( - A coupling model with the good ports to connect to - The binding constraint that ties everything """ + if root.size == 1: + # Nothing to replicate. Just past down the network + return root.network + + master_network = Network(f"{root.id}_{root.network.id}") + + # # We assume that the network was copied already + # # TODO add a guard to verify it - strategy = InvestmentProblemStrategy() - for tree_node in root.traverse(): - ... + # for tree_node in root.traverse(): + # for component_on_node in tree_node.network.components: + # component = component_on_node.replicate(id=f"{tree_node.id}_{component_on_node.id}") + # master_network.add_component(component) - return root.network + return master_network diff --git a/src/andromede/simulation/optimization.py b/src/andromede/simulation/optimization.py index 60cc318b..86a78373 100644 --- a/src/andromede/simulation/optimization.py +++ b/src/andromede/simulation/optimization.py @@ -68,22 +68,7 @@ def _get_parameter_value( ) -> float: data = context.database.get_data(component_id, name) absolute_timestep = context.block_timestep_to_absolute_timestep(block_timestep) - return data.get_value(absolute_timestep, scenario) - - -# TODO: Maybe add the notion of constant parameter in the model -# TODO : And constant over scenarios ? -def _parameter_is_constant_over_time( - component: Component, - name: str, - context: "OptimizationContext", - block_timestep: int, - scenario: int, -) -> bool: - data = context.database.get_data(component.id, name) - return data.get_value(block_timestep, scenario) == IndexingStructure( - time=False, scenario=False - ) + return data.get_value(absolute_timestep, scenario, context.tree_node) class TimestepValueProvider(ABC): @@ -304,12 +289,15 @@ def __init__( block: TimeBlock, scenarios: int, border_management: BlockBorderManagement, + tree_node: str, ): self._network = network self._database = database self._block = block self._scenarios = scenarios self._border_management = border_management + self.tree_node = tree_node + self._component_variables: Dict[TimestepComponentVariableKey, lp.Variable] = {} self._solver_variables: Dict[lp.Variable, SolverVariableInfo] = {} self._connection_fields_expressions: Dict[ @@ -715,6 +703,7 @@ def _register_connection_fields_definitions(self) -> None: ) def _create_variables(self) -> None: + tree_node = self.context.tree_node for component in self.context.network.all_components: component_context = self.context.get_component_context(component) model = component.model @@ -753,7 +742,7 @@ def _create_variables(self) -> None: solver_var = self.solver.NumVar( lower_bound, upper_bound, - f"{component.id}_{model_var.name}_t{block_timestep}_s{scenario}", + f"{tree_node}_{component.id}_{model_var.name}_t{block_timestep}_s{scenario}", ) component_context.add_variable( block_timestep, scenario, model_var.name, solver_var @@ -816,6 +805,7 @@ def build_problem( border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, solver_id: str = "GLOP", problem_strategy: ModelSelectionStrategy = MergedProblemStrategy(), + decision_tree_node: str = "", ) -> OptimizationProblem: """ Entry point to build the optimization problem for a time period. @@ -825,7 +815,7 @@ def build_problem( database.requirements_consistency(network) opt_context = OptimizationContext( - network, database, block, scenarios, border_management + network, database, block, scenarios, border_management, decision_tree_node ) return OptimizationProblem(problem_name, solver, opt_context, problem_strategy) diff --git a/src/andromede/study/data.py b/src/andromede/study/data.py index d1e475f4..db4f0135 100644 --- a/src/andromede/study/data.py +++ b/src/andromede/study/data.py @@ -36,11 +36,7 @@ class ScenarioIndex: @dataclass(frozen=True) class AbstractDataStructure(ABC): @abstractmethod - def get_value( - self, timestep: int, scenario: int, node_id: Optional[int] = None - ) -> ( - float - ): # Is it necessary to add node_id as arguement here ? Yes if TreeData is to be considered as a child class + def get_value(self, timestep: int, scenario: int, node_id: str = "") -> float: """ Get the data value for a given timestep and scenario at a given node Implement this method in subclasses as needed. @@ -60,9 +56,7 @@ def check_requirement(self, time: bool, scenario: bool) -> bool: class ConstantData(AbstractDataStructure): value: float - def get_value( - self, timestep: int, scenario: int, node_id: Optional[int] = None - ) -> float: + def get_value(self, timestep: int, scenario: int, node_id: str = "") -> float: return self.value # ConstantData can be used for time varying or constant models @@ -82,9 +76,7 @@ class TimeSeriesData(AbstractDataStructure): time_series: Mapping[TimeIndex, float] - def get_value( - self, timestep: int, scenario: int, node_id: Optional[int] = None - ) -> float: + def get_value(self, timestep: int, scenario: int, node_id: str = "") -> float: return self.time_series[TimeIndex(timestep)] def check_requirement(self, time: bool, scenario: bool) -> bool: @@ -104,9 +96,7 @@ class ScenarioSeriesData(AbstractDataStructure): scenario_series: Mapping[ScenarioIndex, float] - def get_value( - self, timestep: int, scenario: int, node_id: Optional[int] = None - ) -> float: + def get_value(self, timestep: int, scenario: int, node_id: str = "") -> float: return self.scenario_series[ScenarioIndex(scenario)] def check_requirement(self, time: bool, scenario: bool) -> bool: @@ -126,9 +116,7 @@ class TimeScenarioSeriesData(AbstractDataStructure): time_scenario_series: Mapping[TimeScenarioIndex, float] - def get_value( - self, timestep: int, scenario: int, node_id: Optional[int] = None - ) -> float: + def get_value(self, timestep: int, scenario: int, node_id: str = "") -> float: return self.time_scenario_series[TimeScenarioIndex(timestep, scenario)] def check_requirement(self, time: bool, scenario: bool) -> bool: @@ -140,17 +128,9 @@ def check_requirement(self, time: bool, scenario: bool) -> bool: @dataclass(frozen=True) class TreeData(AbstractDataStructure): - data: Mapping[int, AbstractDataStructure] - - def get_value( - self, timestep: int, scenario: int, node_id: Optional[int] = None - ) -> float: - if ( - not node_id - ): # TODO : Could we remove the default None argument for node_id ? - raise ValueError( - "A node_id must be specified to retrieve a value in TreeData." - ) + data: Mapping[str, AbstractDataStructure] + + def get_value(self, timestep: int, scenario: int, node_id: str = "") -> float: return self.data[node_id].get_value(timestep, scenario) def check_requirement(self, time: bool, scenario: bool) -> bool: @@ -161,7 +141,7 @@ def check_requirement(self, time: bool, scenario: bool) -> bool: @dataclass(frozen=True) -class ComponentParameterIndex: +class DatabaseIndex: component_id: str parameter_name: str @@ -174,22 +154,20 @@ class DataBase: Data can have different structure : constant, varying in time or scenarios. """ - _data: Dict[ComponentParameterIndex, AbstractDataStructure] + _data: Dict[DatabaseIndex, AbstractDataStructure] def __init__(self) -> None: - self._data: Dict[ComponentParameterIndex, AbstractDataStructure] = {} + self._data: Dict[DatabaseIndex, AbstractDataStructure] = {} def get_data(self, component_id: str, parameter_name: str) -> AbstractDataStructure: - return self._data[ComponentParameterIndex(component_id, parameter_name)] + return self._data[DatabaseIndex(component_id, parameter_name)] def add_data( self, component_id: str, parameter_name: str, data: AbstractDataStructure ) -> None: - self._data[ComponentParameterIndex(component_id, parameter_name)] = data + self._data[DatabaseIndex(component_id, parameter_name)] = data - def get_value( - self, index: ComponentParameterIndex, timestep: int, scenario: int - ) -> float: + def get_value(self, index: DatabaseIndex, timestep: int, scenario: int) -> float: if index in self._data: return self._data[index].get_value(timestep, scenario) else: diff --git a/tests/functional/test_investment_pathway.py b/tests/functional/test_investment_pathway.py index 33d4f703..e6ec6674 100644 --- a/tests/functional/test_investment_pathway.py +++ b/tests/functional/test_investment_pathway.py @@ -24,7 +24,7 @@ from andromede.simulation.decision_tree import ( DecisionTreeNode, InterDecisionTimeScenarioConfig, - replicate_network_from_root, + replicate_network_from, ) from andromede.simulation.time_block import TimeBlock from andromede.study.data import ConstantData, DataBase, TreeData @@ -147,13 +147,19 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( # or we index all data, parameters, variables by the resolution node : Make the data struture dependent of the resolution tree... database = DataBase() - database.add_data("N", "spillage", ConstantData(10)) - database.add_data("N", "unsupplied_energy", ConstantData(10000)) + database.add_data("N", "spillage_cost", ConstantData(10)) + database.add_data("N", "ens_cost", ConstantData(10000)) database.add_data( "D", "demand", - TreeData({0: ConstantData(300), 1: ConstantData(600), 2: ConstantData(600)}), + TreeData( + { + "2030": ConstantData(300), + "2040_new_base": ConstantData(600), + "2040_no_base": ConstantData(600), + } + ), ) database.add_data("CAND", "op_cost", ConstantData(10)) @@ -161,18 +167,36 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( database.add_data( "CAND", "invest_cost", - TreeData({0: ConstantData(100), 1: ConstantData(50), 2: ConstantData(50)}), + TreeData( + { + "2030": ConstantData(100), + "2040_new_base": ConstantData(50), + "2040_no_base": ConstantData(50), + } + ), ) database.add_data( "CAND", "max_invest", - TreeData({0: ConstantData(400), 1: ConstantData(100), 2: ConstantData(100)}), + TreeData( + { + "2030": ConstantData(400), + "2040_new_base": ConstantData(100), + "2040_no_base": ConstantData(100), + } + ), ) database.add_data( "BASE", "p_max", - TreeData({0: ConstantData(0), 1: ConstantData(200), 2: ConstantData(0)}), + TreeData( + { + "2030": ConstantData(0), + "2040_new_base": ConstantData(200), + "2040_no_base": ConstantData(0), + } + ), ) database.add_data("BASE", "cost", ConstantData(5)) @@ -198,9 +222,11 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( "2040_no_base", time_scenario_config, parent=decision_tree_root ) - replicate_network_from_root(decision_tree_root) + replicate_network_from(decision_tree_root) decision_coupling_model = model("DECISION_COUPLING") xpansion = build_benders_decomposed_problem( decision_tree_root, database, decision_coupling_model=decision_coupling_model ) + + xpansion.prepare(is_debug=True) diff --git a/tests/integration/test_benders_decomposed.py b/tests/integration/test_benders_decomposed.py index a2ef3269..5c73c55d 100644 --- a/tests/integration/test_benders_decomposed.py +++ b/tests/integration/test_benders_decomposed.py @@ -241,7 +241,7 @@ def test_benders_decomposed_integration( data = { "solution": { "overall_cost": 80_000, - "values": {"CAND_p_max_t0_s0": 100, "DISCRETE_p_max_t0_s0": 100}, + "values": {"_CAND_p_max_t0_s0": 100, "_DISCRETE_p_max_t0_s0": 100}, } } solution = BendersSolution(data) diff --git a/tests/unittests/test_generate_network_on_tree.py b/tests/unittests/test_generate_network_on_tree.py index c77111d3..15ba6088 100644 --- a/tests/unittests/test_generate_network_on_tree.py +++ b/tests/unittests/test_generate_network_on_tree.py @@ -15,7 +15,7 @@ from andromede.simulation.decision_tree import ( DecisionTreeNode, InterDecisionTimeScenarioConfig, - replicate_network_from_root, + replicate_network_from, ) from andromede.study.network import Network @@ -26,19 +26,23 @@ def test_generate_model_on_node() -> None: config = InterDecisionTimeScenarioConfig(blocks, scenarios) network = Network("network_id") - tree_root = DecisionTreeNode("root", config, network) + root = DecisionTreeNode("root", config, network) - assert tree_root.id == "root" - assert tree_root.parent is None - assert not tree_root.children # No children + assert root.id == "root" + assert root.parent is None + assert not root.children # No children - child = DecisionTreeNode("child", config, parent=tree_root) + child = DecisionTreeNode("child", config, parent=root) - print(child.parent) + assert child.parent == root + assert child in root.children - assert child.parent == tree_root - assert child in tree_root.children + grandchild = DecisionTreeNode("grandchild", config, parent=child) - replicate_network_from_root(tree_root) + assert grandchild.parent == child + assert (grandchild not in root.children) and (grandchild in child.children) + + replicate_network_from(root) assert child.network.id == "child_network_id" + assert grandchild.network.id == "grandchild_network_id" From 2171b314b33e1ca3cad58c40b20c098ad31f7003 Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Fri, 26 Apr 2024 10:13:28 +0200 Subject: [PATCH 08/14] Add first pathway test --- src/andromede/libs/standard.py | 10 +- src/andromede/simulation/__init__.py | 1 + .../simulation/benders_decomposed.py | 85 ++++--- src/andromede/simulation/decision_tree.py | 86 +------- src/andromede/simulation/optimization.py | 138 +++++++++--- tests/functional/test_andromede.py | 2 + tests/functional/test_investment_pathway.py | 208 +++++++++++++++--- tests/functional/test_xpansion.py | 2 +- tests/integration/test_benders_decomposed.py | 12 +- .../test_generate_network_on_tree.py | 6 - 10 files changed, 370 insertions(+), 180 deletions(-) diff --git a/src/andromede/libs/standard.py b/src/andromede/libs/standard.py index f9028bf8..48e05d7e 100644 --- a/src/andromede/libs/standard.py +++ b/src/andromede/libs/standard.py @@ -441,10 +441,16 @@ float_variable( "invested_capa", lower_bound=literal(0), - upper_bound=param("max_invest"), structure=CONSTANT, context=ProblemContext.COUPLING, ), + float_variable( + "delta_invest", + lower_bound=literal(0), + upper_bound=param("max_invest"), + structure=CONSTANT, + context=ProblemContext.INVESTMENT, + ), ], ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")], port_fields_definitions=[ @@ -463,5 +469,5 @@ objective_operational_contribution=(param("op_cost") * var("generation")) .sum() .expec(), - objective_investment_contribution=param("invest_cost") * var("invested_capa"), + objective_investment_contribution=param("invest_cost") * var("delta_invest"), ) diff --git a/src/andromede/simulation/__init__.py b/src/andromede/simulation/__init__.py index 23e0a855..46501a29 100644 --- a/src/andromede/simulation/__init__.py +++ b/src/andromede/simulation/__init__.py @@ -14,6 +14,7 @@ BendersDecomposedProblem, build_benders_decomposed_problem, ) +from .decision_tree import DecisionTreeNode, InterDecisionTimeScenarioConfig from .optimization import BlockBorderManagement, OptimizationProblem, build_problem from .output_values import BendersSolution, OutputValues from .runner import BendersRunner, MergeMPSRunner diff --git a/src/andromede/simulation/benders_decomposed.py b/src/andromede/simulation/benders_decomposed.py index 5e42306f..842dd10b 100644 --- a/src/andromede/simulation/benders_decomposed.py +++ b/src/andromede/simulation/benders_decomposed.py @@ -18,14 +18,12 @@ import pathlib from typing import Any, Dict, List, Optional -from anytree import Node as TreeNode - -from andromede.model.model import Model -from andromede.simulation.decision_tree import DecisionTreeNode, create_master_network +from andromede.simulation.decision_tree import DecisionTreeNode from andromede.simulation.optimization import ( BlockBorderManagement, OptimizationProblem, build_problem, + fusion_problems, ) from andromede.simulation.output_values import ( BendersDecomposedSolution, @@ -79,7 +77,6 @@ def export_structure(self) -> str: """ 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 @@ -89,11 +86,10 @@ def export_structure(self) -> str: 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) + 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] = {} @@ -126,7 +122,7 @@ def export_options( "TRACE": True, "SLAVE_WEIGHT": "CONSTANT", "SLAVE_WEIGHT_VALUE": 1, - "MASTER_NAME": "master", + "MASTER_NAME": f"{self.master.name}", "STRUCTURE_FILE": "structure.txt", "INPUTROOT": ".", "CSV_NAME": "benders_output_trace", @@ -144,19 +140,21 @@ def export_options( } return options_values - def prepare( + def initialise( self, *, solver_name: str = "XPRESS", log_level: int = 0, is_debug: bool = False, ) -> None: - serialize("master.mps", self.master.export_as_mps(), self.emplacement) + serialize( + f"{self.master.name}.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(f"structure.txt", self.export_structure(), self.emplacement) serialize_json( "options.json", self.export_options(solver_name=solver_name, log_level=log_level), @@ -164,7 +162,9 @@ def prepare( ) if is_debug: - serialize("master.lp", self.master.export_as_lp(), self.emplacement) + serialize( + f"{self.master.name}.lp", self.master.export_as_lp(), self.emplacement + ) for subproblem in self.subproblems: serialize( f"{subproblem.name}.lp", subproblem.export_as_lp(), self.emplacement @@ -192,7 +192,7 @@ def run( log_level: int = 0, should_merge: bool = False, ) -> bool: - self.prepare(solver_name=solver_name, log_level=log_level) + self.initialise(solver_name=solver_name, log_level=log_level) if not should_merge: return_code = BendersRunner(self.emplacement).run() @@ -211,9 +211,9 @@ def build_benders_decomposed_problem( decision_tree_root: DecisionTreeNode, database: DataBase, *, - decision_coupling_model: Optional[Model] = None, border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, solver_id: str = "GLOP", + coupling_network: Network = Network(""), ) -> BendersDecomposedProblem: """ Entry point to build the xpansion problem for a time period @@ -221,38 +221,57 @@ def build_benders_decomposed_problem( Returns a Benders Decomposed problem """ - master_network = create_master_network(decision_tree_root, decision_coupling_model) + null_time_block = TimeBlock( + 0, [0] + ) # Not necessary for master, but list must be non-empty + null_scenario = 0 # Not necessary for master - # Benders Decomposed Master Problem - master = build_problem( - master_network, + coupler = build_problem( + coupling_network, database, - 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, + null_time_block, + null_scenario, + problem_name="coupler", solver_id=solver_id, - problem_strategy=InvestmentProblemStrategy(), + build_strategy=InvestmentProblemStrategy(), ) - # Benders Decomposed Sub-problems - subproblems = [] + masters = [] # Benders Decomposed Master Problem + subproblems = [] # Benders Decomposed Sub-problems + for tree_node in decision_tree_root.traverse(): + suffix = f"_{tree_node.id}" if decision_tree_root.size > 1 else "" + + masters.append( + build_problem( + tree_node.network, + database, + null_time_block, + null_scenario, + problem_name=f"master{suffix}", + solver_id=solver_id, + build_strategy=InvestmentProblemStrategy(), + decision_tree_node=tree_node.id, + ) + ) + for block in tree_node.config.blocks: - # Xpansion Sub-problems + if len(tree_node.config.blocks) > 1: + suffix += f"_t{block.id}" + subproblems.append( build_problem( tree_node.network, database, block, tree_node.config.scenarios, - problem_name=f"subproblem_{tree_node.id}_{block.id}", + problem_name=f"subproblem{suffix}", solver_id=solver_id, - problem_strategy=OperationalProblemStrategy(), + build_strategy=OperationalProblemStrategy(), decision_tree_node=tree_node.id, ) ) + master = fusion_problems(masters, coupler) + return BendersDecomposedProblem(master, subproblems) diff --git a/src/andromede/simulation/decision_tree.py b/src/andromede/simulation/decision_tree.py index 505a6904..c078f1d7 100644 --- a/src/andromede/simulation/decision_tree.py +++ b/src/andromede/simulation/decision_tree.py @@ -15,8 +15,6 @@ from anytree import LevelOrderIter, NodeMixin -from andromede.model.model import Model -from andromede.simulation.strategy import InvestmentProblemStrategy from andromede.simulation.time_block import TimeBlock from andromede.study.network import Network @@ -47,83 +45,7 @@ def __init__( if children: self.children = children - def traverse(self) -> Generator["DecisionTreeNode", None, None]: - yield from LevelOrderIter(self) - - -def replicate_network_from(root: DecisionTreeNode) -> None: - if root.size == 1: - # Nothing to replicate. Just past down the network - return - - else: - # Replicates the network and updates their id to take into account the tree node's id - original_network_id = root.network.id - for tree_node in root.traverse(): - tree_node.network = root.network.replicate( - id=f"{tree_node.id}_{original_network_id}" - ) - - -def create_master_network( - root: DecisionTreeNode, - decision_coupling_model: Optional[Model], -) -> Network: - # TODO Use ports for coupling different models across the decision tree - """ - Each candidate model should have one of these ports. - As for balance, they all should have a common name like "investment" and they could be defined as: - - ports=[ModelPort(port_type=INVESTMENT_PORT_TYPE, port_name="investment_port")], - port_fields_definitions=[ - PortFieldDefinition( - port_field=PortFieldId("investment_port", "investment"), - definition=param("p_max"), ---> p_max if we are investing the p_max for instance - ) - ], - - Then, here, we would traverse the tree and create a Network that only keeps - variables allowed in the InvestmentProblemStrategy (we would need to change their names to something like - p_max_root, p_max_child_id, p_max_child_id2, ...) - We would connect them by some constraints using these ports to the decision_coupling_model (a parameter here), - so we could have something like: - - master_network.connect(PortRef(candidate_on_root, "investment_port"), PortRef(coupling_model, "investment_port")) - master_network.connect(PortRef(candidate_on_child, "investment_port"), PortRef(coupling_model, "investment_port")) - - On the coupling model, we would have something like for the nodes: - ports=[ModelPort(port_type=INVESTMENT_PORT_TYPE, port_name="investment_port")], - binding_constraints=[ - Constraint( - name="Pathway_Investment", - expression=port_field("investment_port", "flow").some_operator(), - ) - ], - - Maybe we would have to define a operator to represent consecutive inequalities ? - Or create one coupling model per pair parent-child so the expression would become - port_field("investment_port", "flow").sum_connections() <= 0 ? In this case, we would have to define - a negative and a positive value - - To resume, a network here would need: - - All candidates on all tree nodes (so we don't make great changes on build_problem) - - Update investment variable names to show their corresponding tree node; - - Add a investment port to each model - - A coupling model with the good ports to connect to - - The binding constraint that ties everything - """ - if root.size == 1: - # Nothing to replicate. Just past down the network - return root.network - - master_network = Network(f"{root.id}_{root.network.id}") - - # # We assume that the network was copied already - # # TODO add a guard to verify it - - # for tree_node in root.traverse(): - # for component_on_node in tree_node.network.components: - # component = component_on_node.replicate(id=f"{tree_node.id}_{component_on_node.id}") - # master_network.add_component(component) - - return master_network + def traverse( + self, depth: Optional[int] = None + ) -> Generator["DecisionTreeNode", None, None]: + yield from LevelOrderIter(self, maxlevel=depth) diff --git a/src/andromede/simulation/optimization.py b/src/andromede/simulation/optimization.py index 86a78373..bef4ad54 100644 --- a/src/andromede/simulation/optimization.py +++ b/src/andromede/simulation/optimization.py @@ -289,17 +289,19 @@ def __init__( block: TimeBlock, scenarios: int, border_management: BlockBorderManagement, - tree_node: str, + build_strategy: ModelSelectionStrategy = MergedProblemStrategy(), + decision_tree_node: str = "", ): self._network = network self._database = database self._block = block self._scenarios = scenarios self._border_management = border_management - self.tree_node = tree_node + self._build_strategy = build_strategy + self._tree_node = decision_tree_node self._component_variables: Dict[TimestepComponentVariableKey, lp.Variable] = {} - self._solver_variables: Dict[lp.Variable, SolverVariableInfo] = {} + self._solver_variables: Dict[str, SolverVariableInfo] = {} self._connection_fields_expressions: Dict[ PortFieldKey, List[ExpressionNode] ] = {} @@ -312,6 +314,14 @@ def network(self) -> Network: def scenarios(self) -> int: return self._scenarios + @property + def tree_node(self) -> str: + return self._tree_node + + @property + def strategy(self) -> ModelSelectionStrategy: + return self._build_strategy + def block_length(self) -> int: return len(self._block.timesteps) @@ -372,19 +382,19 @@ def register_component_variable( block_timestep: int, scenario: int, component_id: str, - variable_name: str, + model_var_name: str, variable: lp.Variable, ) -> None: key = TimestepComponentVariableKey( - component_id, variable_name, block_timestep, scenario + component_id, model_var_name, block_timestep, scenario ) if key not in self._component_variables: - self._solver_variables[variable] = SolverVariableInfo( + self._solver_variables[variable.name()] = SolverVariableInfo( variable.name(), len(self._solver_variables), False ) self._component_variables[key] = variable - def get_component_context(self, component: Component) -> ComponentContext: + def create_component_context(self, component: Component) -> ComponentContext: return ComponentContext(self, component) def register_connection_fields_expressions( @@ -511,7 +521,7 @@ def _create_objective( ) for solver_var in solver_vars: - opt_context._solver_variables[solver_var].is_in_objective = True + opt_context._solver_variables[solver_var.name()].is_in_objective = True obj.SetCoefficient( solver_var, obj.GetCoefficient(solver_var) + weight * term.coefficient, @@ -655,19 +665,16 @@ class OptimizationProblem: name: str solver: lp.Solver context: OptimizationContext - strategy: ModelSelectionStrategy 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() @@ -703,12 +710,11 @@ def _register_connection_fields_definitions(self) -> None: ) def _create_variables(self) -> None: - tree_node = self.context.tree_node for component in self.context.network.all_components: - component_context = self.context.get_component_context(component) + component_context = self.context.create_component_context(component) model = component.model - for model_var in self.strategy.get_variables(model): + for model_var in self.context.strategy.get_variables(model): var_indexing = IndexingStructure( model_var.structure.time, model_var.structure.scenario ) @@ -722,7 +728,21 @@ def _create_variables(self) -> None: instantiated_ub_expr = _instantiate_model_expression( model_var.upper_bound, component.id, self.context ) + + # Set solver var name + # 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_name: str = f"{model_var.name}" + if component.id: + solver_var_name = f"{component.id}_" + solver_var_name + if self.context.tree_node: + solver_var_name = f"{self.context.tree_node}_" + solver_var_name + for block_timestep in self.context.get_time_indices(var_indexing): + if self.context.block_length() > 1: + solver_var_name = solver_var_name + f"_t{block_timestep}" + for scenario in self.context.get_scenario_indices(var_indexing): lower_bound = -self.solver.infinity() upper_bound = self.solver.infinity() @@ -735,14 +755,14 @@ def _create_variables(self) -> None: instantiated_ub_expr ).get_value(block_timestep, scenario) + if self.context.scenarios > 1: + solver_var_name = solver_var_name + f"_s{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"{tree_node}_{component.id}_{model_var.name}_t{block_timestep}_s{scenario}", + solver_var_name, ) component_context.add_variable( block_timestep, scenario, model_var.name, solver_var @@ -750,7 +770,7 @@ def _create_variables(self) -> None: def _create_constraints(self) -> None: for component in self.context.network.all_components: - for constraint in self.strategy.get_constraints(component.model): + for constraint in self.context.strategy.get_constraints(component.model): instantiated_expr = _instantiate_model_expression( constraint.expression, component.id, self.context ) @@ -769,16 +789,16 @@ def _create_constraints(self) -> None: ) _create_constraint( self.solver, - self.context.get_component_context(component), + self.context.create_component_context(component), instantiated_constraint, ) def _create_objectives(self) -> None: for component in self.context.network.all_components: - component_context = self.context.get_component_context(component) + component_context = self.context.create_component_context(component) model = component.model - for objective in self.strategy.get_objectives(model): + for objective in self.context.strategy.get_objectives(model): if objective is not None: _create_objective( self.solver, @@ -804,7 +824,7 @@ def build_problem( problem_name: str = "optimization_problem", border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, solver_id: str = "GLOP", - problem_strategy: ModelSelectionStrategy = MergedProblemStrategy(), + build_strategy: ModelSelectionStrategy = MergedProblemStrategy(), decision_tree_node: str = "", ) -> OptimizationProblem: """ @@ -815,7 +835,75 @@ def build_problem( database.requirements_consistency(network) opt_context = OptimizationContext( - network, database, block, scenarios, border_management, decision_tree_node + network, + database, + block, + scenarios, + border_management, + build_strategy, + decision_tree_node, ) - return OptimizationProblem(problem_name, solver, opt_context, problem_strategy) + return OptimizationProblem(problem_name, solver, opt_context) + + +def fusion_problems( + masters: List[OptimizationProblem], coupler: OptimizationProblem +) -> OptimizationProblem: + if len(masters) == 1: + # Nothing to fusion. Just past down the master + return masters[0] + + root_master = coupler + root_master.name = "master" + + root_vars: Dict[str, lp.Variable] = dict() + root_cstrs: Dict[str, lp.Constraint] = dict() + root_objective = root_master.solver.Objective() + + # We stock the coupler's variables to check for + # same name variables in the masters + for var in root_master.solver.variables(): + root_vars[var.name()] = var + + for master in masters: + context = master.context + objective = master.solver.Objective() + + for var in master.solver.variables(): + # If variable not already in coupler, we add it + # Otherwise we update its upper and lower bounds + if var.name() not in root_vars: + root_var = root_master.solver.NumVar(var.lb(), var.ub(), var.name()) + root_master.context._solver_variables[var.name()] = SolverVariableInfo( + var.name(), + len(root_master.context._solver_variables), + context._solver_variables[var.name()].is_in_objective, + ) + else: + root_var = root_vars[var.name()] + root_var.SetLb(var.lb()) + root_var.SetUb(var.ub()) + root_master.context._solver_variables[ + var.name() + ].is_in_objective = context._solver_variables[ + var.name() + ].is_in_objective + + for cstr in master.solver.constraints(): + coeff = cstr.GetCoefficient(var) + # If variable present in constraint, we add the constraint to root + if coeff != 0: + key = f"{master.name}_{cstr.name()}" + if key not in root_cstrs: + root_cstrs[key] = root_master.solver.Constraint( + cstr.Lb(), cstr.Ub(), key + ) + root_cstr = root_cstrs[key] + root_cstr.SetCoefficient(root_var, coeff) + + obj_coeff = objective.GetCoefficient(var) + if obj_coeff != 0: + root_objective.SetCoefficient(root_var, obj_coeff) + + return root_master diff --git a/tests/functional/test_andromede.py b/tests/functional/test_andromede.py index 300ff88d..ba033a32 100644 --- a/tests/functional/test_andromede.py +++ b/tests/functional/test_andromede.py @@ -30,6 +30,8 @@ from andromede.model.model import PortFieldDefinition, PortFieldId from andromede.simulation import ( BlockBorderManagement, + DecisionTreeNode, + InterDecisionTimeScenarioConfig, OutputValues, TimeBlock, build_problem, diff --git a/tests/functional/test_investment_pathway.py b/tests/functional/test_investment_pathway.py index e6ec6674..a1599ec3 100644 --- a/tests/functional/test_investment_pathway.py +++ b/tests/functional/test_investment_pathway.py @@ -13,20 +13,27 @@ import pytest from anytree import Node as TreeNode +from andromede.expression import literal, var +from andromede.expression.indexing_structure import IndexingStructure from andromede.libs.standard import ( DEMAND_MODEL, GENERATOR_MODEL, NODE_WITH_SPILL_AND_ENS, THERMAL_CANDIDATE_WITH_ALREADY_INSTALLED_CAPA, ) +from andromede.model.common import ProblemContext +from andromede.model.constraint import Constraint from andromede.model.model import model -from andromede.simulation.benders_decomposed import build_benders_decomposed_problem +from andromede.model.variable import float_variable +from andromede.simulation import ( + BendersSolution, + TimeBlock, + build_benders_decomposed_problem, +) from andromede.simulation.decision_tree import ( DecisionTreeNode, InterDecisionTimeScenarioConfig, - replicate_network_from, ) -from andromede.simulation.time_block import TimeBlock from andromede.study.data import ConstantData, DataBase, TreeData from andromede.study.network import Component, Network, Node, PortRef, create_component @@ -60,6 +67,158 @@ def node() -> Node: return node +def test_investment_pathway_on_sequential_nodes( + node: Node, + demand: Component, + candidate: Component, +) -> None: + database = DataBase() + database.add_data("N", "spillage_cost", ConstantData(10)) + database.add_data("N", "ens_cost", ConstantData(10000)) + + database.add_data( + "D", + "demand", + TreeData( + { + "parent": ConstantData(100), + "child": ConstantData(200), + } + ), + ) + + database.add_data("CAND", "op_cost", ConstantData(10)) + database.add_data("CAND", "already_installed_capa", ConstantData(100)) + + database.add_data( + "CAND", + "invest_cost", + TreeData( + { + "parent": ConstantData(100), + "child": ConstantData(300), + } + ), + ) + + database.add_data( + "CAND", + "max_invest", + TreeData( + { + "parent": ConstantData(80), + "child": ConstantData(100), + } + ), + ) + + COUPLING_MODEL = model( + id="COUPLING", + variables=[ + float_variable( + "parent_CAND_delta_invest", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + float_variable( + "parent_CAND_invested_capa", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + float_variable( + "child_CAND_delta_invest", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + float_variable( + "child_CAND_invested_capa", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + ], + constraints=[ + Constraint( + name="Max investment on parent", + expression=var("parent_CAND_invested_capa") + == var("parent_CAND_delta_invest"), + context=ProblemContext.INVESTMENT, + ), + Constraint( + name="Max investment on child", + expression=var("child_CAND_invested_capa") + == var("child_CAND_delta_invest") + var("parent_CAND_invested_capa"), + context=ProblemContext.INVESTMENT, + ), + ], + ) + + network_coupling = Network("coupling_test") + network_coupling.add_component(create_component(model=COUPLING_MODEL, id="")) + + demand_par = demand.replicate() + candidate_par = candidate.replicate() + + network_par = Network("parent_test") + network_par.add_node(node) + network_par.add_component(demand_par) + network_par.add_component(candidate_par) + network_par.connect( + PortRef(demand_par, "balance_port"), PortRef(node, "balance_port") + ) + network_par.connect( + PortRef(candidate_par, "balance_port"), PortRef(node, "balance_port") + ) + + demand_chd = demand.replicate() + candidate_chd = candidate.replicate() + + network_chd = Network("child_test") + network_chd.add_node(node) + network_chd.add_component(demand_chd) + network_chd.add_component(candidate_chd) + network_chd.connect( + PortRef(demand_chd, "balance_port"), PortRef(node, "balance_port") + ) + network_chd.connect( + PortRef(candidate_chd, "balance_port"), PortRef(node, "balance_port") + ) + + config = InterDecisionTimeScenarioConfig([TimeBlock(0, [0])], 1) + + decision_tree_par = DecisionTreeNode("parent", config, network_par) + decision_tree_chd = DecisionTreeNode( + "child", config, network_chd, parent=decision_tree_par + ) + + xpansion = build_benders_decomposed_problem( + decision_tree_par, database, coupling_network=network_coupling + ) + + data = { + "solution": { + "overall_cost": 17_000, + "values": { + "parent_CAND_delta_invest": 80, + "child_CAND_delta_invest": 20, + "parent_CAND_invested_capa": 80, + "child_CAND_invested_capa": 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}" + + def test_investment_pathway_on_a_tree_with_one_root_two_children( generator: Component, candidate: Component, @@ -155,9 +314,9 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( "demand", TreeData( { - "2030": ConstantData(300), - "2040_new_base": ConstantData(600), - "2040_no_base": ConstantData(600), + "ROOT": ConstantData(300), + "CHILD_A": ConstantData(600), + "CHILD_B": ConstantData(600), } ), ) @@ -169,9 +328,9 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( "invest_cost", TreeData( { - "2030": ConstantData(100), - "2040_new_base": ConstantData(50), - "2040_no_base": ConstantData(50), + "ROOT": ConstantData(100), + "CHILD_A": ConstantData(50), + "CHILD_B": ConstantData(50), } ), ) @@ -180,9 +339,9 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( "max_invest", TreeData( { - "2030": ConstantData(400), - "2040_new_base": ConstantData(100), - "2040_no_base": ConstantData(100), + "ROOT": ConstantData(400), + "CHILD_A": ConstantData(100), + "CHILD_B": ConstantData(100), } ), ) @@ -192,9 +351,9 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( "p_max", TreeData( { - "2030": ConstantData(0), - "2040_new_base": ConstantData(200), - "2040_no_base": ConstantData(0), + "ROOT": ConstantData(0), + "CHILD_A": ConstantData(200), + "CHILD_B": ConstantData(0), } ), ) @@ -214,19 +373,14 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( [TimeBlock(0, [0])], scenarios ) - decision_tree_root = DecisionTreeNode("2030", time_scenario_config, network) - decision_tree_new_base = DecisionTreeNode( - "2040_new_base", time_scenario_config, parent=decision_tree_root + dt_root = DecisionTreeNode("ROOT", time_scenario_config, network) + dt_child_A = DecisionTreeNode( + "CHILD_A", time_scenario_config, network, parent=dt_root ) - decision_tree_no_base = DecisionTreeNode( - "2040_no_base", time_scenario_config, parent=decision_tree_root + dt_child_B = DecisionTreeNode( + "CHILD_B", time_scenario_config, network, parent=dt_root ) - replicate_network_from(decision_tree_root) - decision_coupling_model = model("DECISION_COUPLING") - - xpansion = build_benders_decomposed_problem( - decision_tree_root, database, decision_coupling_model=decision_coupling_model - ) + xpansion = build_benders_decomposed_problem(dt_root, database) - xpansion.prepare(is_debug=True) + # xpansion.initialise(is_debug=True) diff --git a/tests/functional/test_xpansion.py b/tests/functional/test_xpansion.py index 364a1a05..b36307b3 100644 --- a/tests/functional/test_xpansion.py +++ b/tests/functional/test_xpansion.py @@ -185,7 +185,7 @@ def test_generation_xpansion_single_time_step_single_scenario( database, TimeBlock(1, [0]), scenarios, - problem_strategy=MergedProblemStrategy(), + build_strategy=MergedProblemStrategy(), ) status = problem.solver.Solve() diff --git a/tests/integration/test_benders_decomposed.py b/tests/integration/test_benders_decomposed.py index 5c73c55d..1990029a 100644 --- a/tests/integration/test_benders_decomposed.py +++ b/tests/integration/test_benders_decomposed.py @@ -234,14 +234,18 @@ def test_benders_decomposed_integration( blocks = [TimeBlock(1, [0])] config = InterDecisionTimeScenarioConfig(blocks, scenarios) - decision_tree_root = DecisionTreeNode("root", config, network) + decision_tree_root = DecisionTreeNode("", config, network) xpansion = build_benders_decomposed_problem(decision_tree_root, database) data = { "solution": { "overall_cost": 80_000, - "values": {"_CAND_p_max_t0_s0": 100, "_DISCRETE_p_max_t0_s0": 100}, + "values": { + "CAND_p_max": 100, + "DISCRETE_p_max": 100, + "DISCRETE_nb_units": 10, + }, } } solution = BendersSolution(data) @@ -324,7 +328,7 @@ def test_benders_decomposed_multi_time_block_single_scenario( blocks = [TimeBlock(1, [0]), TimeBlock(2, [1])] config = InterDecisionTimeScenarioConfig(blocks, scenarios) - decision_tree_root = DecisionTreeNode("root", config, network) + decision_tree_root = DecisionTreeNode("", config, network) xpansion = build_benders_decomposed_problem(decision_tree_root, database) @@ -332,7 +336,7 @@ def test_benders_decomposed_multi_time_block_single_scenario( "solution": { "overall_cost": 62_000, "values": { - "CAND_p_max_t0_s0": 100, + "CAND_p_max": 100, }, } } diff --git a/tests/unittests/test_generate_network_on_tree.py b/tests/unittests/test_generate_network_on_tree.py index 15ba6088..20d93a0b 100644 --- a/tests/unittests/test_generate_network_on_tree.py +++ b/tests/unittests/test_generate_network_on_tree.py @@ -15,7 +15,6 @@ from andromede.simulation.decision_tree import ( DecisionTreeNode, InterDecisionTimeScenarioConfig, - replicate_network_from, ) from andromede.study.network import Network @@ -41,8 +40,3 @@ def test_generate_model_on_node() -> None: assert grandchild.parent == child assert (grandchild not in root.children) and (grandchild in child.children) - - replicate_network_from(root) - - assert child.network.id == "child_network_id" - assert grandchild.network.id == "grandchild_network_id" From a76b96b3eed44ac436297d56e1e55ba8b38dcaea Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Thu, 30 May 2024 17:22:24 +0200 Subject: [PATCH 09/14] Add three node test for pathway --- tests/functional/test_investment_pathway.py | 219 ++++++++++++++------ 1 file changed, 156 insertions(+), 63 deletions(-) diff --git a/tests/functional/test_investment_pathway.py b/tests/functional/test_investment_pathway.py index a1599ec3..5774cae7 100644 --- a/tests/functional/test_investment_pathway.py +++ b/tests/functional/test_investment_pathway.py @@ -226,9 +226,12 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( node: Node, ) -> None: """ - This use case aims at representing the situation where investment decisions are to be made at different, say "planning times". An actualisation rate can be taken into account. + This use case aims at representing the situation where investment decisions are to be made at different, say "planning times". + An actualization rate can be taken into account. - The novelty compared the actual usage of planning tools, is that the planning decisions at a given time are taken without knowing exactly which "macro-scenario" / hypothesis on the system that will eventually happen (only knowing the probability distribution of these hypothesis). + The novelty compared the actual usage of planning tools, is that the planning decisions at a given time + are taken without knowing exactly which "macro-scenario" / hypothesis on the system that will eventually happen + (only knowing the probability distribution of these hypothesis). This example models a case where investment decisions have to be made in 2030 and 2040. - In 2030, we have full knowledge of the existing assets @@ -236,10 +239,10 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( - A case where there is no change in the generation assets since 2030 (except te potential investment in 2030) - A case where a base generation unit is present - When taking the decision in 2030, we do not know which case will occur in 2040 and we seek the best decision given a risk criterion (the expectation here). - - The value of these models lies in the output for the first decision rather than the decisions at the later stages as the first decisions are related to "what we have to do today" ? - + When taking the decision in 2030, we do not know which case will occur in 2040 + and we seek the best decision given a risk criterion (the expectation here). + The value of these models lies in the output for the first decision rather than + the decisions at the later stages as the first decisions are related to "what we have to do today" ? More specifically, to define the use case, we define the following tree representing the system at the different decision times and hypothesis 2030 (root node) : @@ -250,9 +253,9 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( Max investment = 400, Investment cost = 100 Unsupplied energy : - Cost = 10000 + Cost = 10 000 - 2040 with new base (scenario 1) : + 2040 with new base (child A) : Demand = 600 Generator : P_max = 200, @@ -263,9 +266,9 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( P_max = 200, Production cost = 5 Unsupplied energy : - Cost = 10000 + Cost = 10 000 - 2040 no base (scenario 2) : + 2040 no base (child B) : Demand = 600 Generator : P_max = 200, @@ -273,50 +276,54 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( Max investment = 100, Investment cost = 50 Unsupplied energy : - Cost = 10000 + Cost = 10 000 - In the second decision time, demand increases from 300 to 600 in both scenarios. However, investment capacity in the candidate is limited to 100 in the second stage. Investment cost decreases to reflect the effect of a discount rate. + In the second decision time, demand increases from 300 to 600 in both scenarios. + However, investment capacity in the candidate is limited to 100 in the second stage. + Investment cost decreases to reflect the effect of a discount rate. - In case 1, a base unit of capacity 100 has arrived and can produce at the same cost than the candidate. As it is more intersting to invest the latest possible, the optimal solution for scenario 1 is to invest [100, 100]. + In case 1, a base unit of capacity 100 has arrived and can produce at smaller cost than the candidate. + As it is more interesting to invest the latest possible, the optimal solution for this scenario is to invest [100, 100]. - In case 2, there is no base unit and the max investment is 100 in the second stage, therefore if we consider scenario 2 only, as unsupplied energy is very expensive, the best investment is [300, 100] + In case 2, there is no base unit and the max investment is 100 in the second stage, + therefore if we consider scenario 2 only, as unsupplied energy is very expensive, the best investment is [300, 100] But here as we solve on the tree, we need to find the best solution in expectation on the set of paths in the tree. - With initial investment = 100 : - Total cost = [100 x 100 (investment root) + 10 x 300 (prod root)] - + 0.5 (proba child 1) x [100 x 50 (investment child 1) + 10 x 400 (prod generator) + 5 x 200 (prod base)] - + 0.5 (proba child 2) x [100 x 50 (investment child 2) + 10 x 400 (prod generator) + 1000 x 200 (unsupplied energy)] - = 122 500 + Case 1 : prob | investment | operational + root: 1 x [ 100 x 100 + 10 x 300 ] + child A: + 0.5 x [ 50 x 100 + 10 x 400 (generator) + 5 x 200 (base)] + child B: + 0.5 x [ 50 x 100 + 10 x 400 (generator) + 10 000 x 200 (unsupplied energy)] + = 1 022 500 - With initial investment = 300 : - Total cost = [100 x 300 (investment root) + 10 x 300 (prod root)] - + 0.5 (proba child 1) x [10 x 400 (prod generator) + 5 x 200 (prod base)] - + 0.5 (proba child 2) x [100 x 50 (investment child 2) + 10 x 600 (prod generator)] - = 41 000 + Case 2 : prob | investment | operational + root: 1 x [ 100 x 300 + 10 x 300 ] + child A: + 0.5 x [ 50 x 0 + 10 x 400 (generator) + 5 x 200 (base)] + child B: + 0.5 x [ 50 x 100 + 10 x 600 (generator)] + = 41 000 - As investing less than 300 in the first stage would increase the unsupplied energy and lead to an increase in overall cost (-1 MW invested in 1st stage => + 1 MW unsp energy => +900/MW cost increase more or less), the optimal solution is to invest : + As investing less than 300 in the first stage would increase the unsupplied energy and lead to an increase in overall cost + (-1 MW invested in 1st stage => + 1 MW unsupplied energy => +900/MW cost increase more or less), the optimal solution is to invest : - 300 at first stage - - 0 in child 1 - - 100 in child 2 - + - 0 in child A + - 100 in child B """ # Either we duplicate all network for each node : Lots of duplications - # or we index all data, parameters, variables by the resolution node : Make the data struture dependent of the resolution tree... + # or we index all data, parameters, variables by the resolution node : Make the data structure dependent of the resolution tree... database = DataBase() database.add_data("N", "spillage_cost", ConstantData(10)) - database.add_data("N", "ens_cost", ConstantData(10000)) + database.add_data("N", "ens_cost", ConstantData(10_000)) database.add_data( "D", "demand", TreeData( { - "ROOT": ConstantData(300), - "CHILD_A": ConstantData(600), - "CHILD_B": ConstantData(600), + "root": ConstantData(300), + "childA": ConstantData(600), + "childB": ConstantData(600), } ), ) @@ -328,9 +335,9 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( "invest_cost", TreeData( { - "ROOT": ConstantData(100), - "CHILD_A": ConstantData(50), - "CHILD_B": ConstantData(50), + "root": ConstantData(100), + "childA": ConstantData(50), + "childB": ConstantData(50), } ), ) @@ -339,48 +346,134 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( "max_invest", TreeData( { - "ROOT": ConstantData(400), - "CHILD_A": ConstantData(100), - "CHILD_B": ConstantData(100), + "root": ConstantData(400), + "childA": ConstantData(100), + "childB": ConstantData(100), } ), ) - database.add_data( - "BASE", - "p_max", - TreeData( - { - "ROOT": ConstantData(0), - "CHILD_A": ConstantData(200), - "CHILD_B": ConstantData(0), - } - ), - ) + database.add_data("BASE", "p_max", ConstantData(200)) database.add_data("BASE", "cost", ConstantData(5)) - 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")) + COUPLING_MODEL = model( + id="COUPLING", + variables=[ + float_variable( + "root_CAND_delta_invest", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + float_variable( + "root_CAND_invested_capa", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + float_variable( + "childA_CAND_delta_invest", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + float_variable( + "childA_CAND_invested_capa", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + float_variable( + "childB_CAND_delta_invest", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + float_variable( + "childB_CAND_invested_capa", + lower_bound=literal(0), + structure=IndexingStructure(False, False), + context=ProblemContext.INVESTMENT, + ), + ], + constraints=[ + Constraint( + name="Max investment on root", + expression=var("root_CAND_invested_capa") + == var("root_CAND_delta_invest"), + context=ProblemContext.INVESTMENT, + ), + Constraint( + name="Max investment on child A", + expression=var("childA_CAND_invested_capa") + == var("childA_CAND_delta_invest") + var("root_CAND_invested_capa"), + context=ProblemContext.INVESTMENT, + ), + Constraint( + name="Max investment on child B", + expression=var("childB_CAND_invested_capa") + == var("childB_CAND_delta_invest") + var("root_CAND_invested_capa"), + context=ProblemContext.INVESTMENT, + ), + ], + ) + + network_coupler = Network("coupling_test") + network_coupler.add_component(create_component(model=COUPLING_MODEL, id="")) + + network_root = Network("root_network") + network_root.add_node(node) + network_root.add_component(demand) + network_root.add_component(candidate) + network_root.connect(PortRef(demand, "balance_port"), PortRef(node, "balance_port")) + network_root.connect( + PortRef(candidate, "balance_port"), PortRef(node, "balance_port") + ) + + network_childA = network_root.replicate(id="childA_network") + network_childA.add_component(generator) + network_childA.connect( + PortRef(generator, "balance_port"), PortRef(node, "balance_port") + ) + + network_childB = network_root.replicate(id="childB_network") scenarios = 1 time_scenario_config = InterDecisionTimeScenarioConfig( [TimeBlock(0, [0])], scenarios ) - dt_root = DecisionTreeNode("ROOT", time_scenario_config, network) + dt_root = DecisionTreeNode("root", time_scenario_config, network_root) dt_child_A = DecisionTreeNode( - "CHILD_A", time_scenario_config, network, parent=dt_root + "childA", time_scenario_config, network_childA, parent=dt_root ) dt_child_B = DecisionTreeNode( - "CHILD_B", time_scenario_config, network, parent=dt_root + "childB", time_scenario_config, network_childB, parent=dt_root ) - xpansion = build_benders_decomposed_problem(dt_root, database) + xpansion = build_benders_decomposed_problem( + dt_root, database, coupling_network=network_coupler + ) + xpansion.initialise(is_debug=True) - # xpansion.initialise(is_debug=True) + data = { + "solution": { + "overall_cost": 49_000, + "values": { + "root_CAND_delta_invest": 300, + "childA_CAND_delta_invest": 0, + "childB_CAND_delta_invest": 100, + "root_CAND_invested_capa": 300, + "childA_CAND_invested_capa": 300, + "childB_CAND_invested_capa": 400, + }, + } + } + 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}" From 1e51b36b8900b77288b1462e9d391d1f99125fe9 Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Fri, 31 May 2024 17:00:01 +0200 Subject: [PATCH 10/14] Add prob to tree nodes for expectation computation --- src/andromede/simulation/decision_tree.py | 7 ++++++ tests/functional/test_investment_pathway.py | 23 ++++++++++--------- .../test_generate_network_on_tree.py | 15 +++++++++++- 3 files changed, 33 insertions(+), 12 deletions(-) diff --git a/src/andromede/simulation/decision_tree.py b/src/andromede/simulation/decision_tree.py index c078f1d7..f9404e5d 100644 --- a/src/andromede/simulation/decision_tree.py +++ b/src/andromede/simulation/decision_tree.py @@ -29,6 +29,7 @@ class DecisionTreeNode(NodeMixin): id: str config: InterDecisionTimeScenarioConfig network: Network + prob: float def __init__( self, @@ -37,11 +38,17 @@ def __init__( network: Network = Network(""), parent: Optional["DecisionTreeNode"] = None, children: Optional[Iterable["DecisionTreeNode"]] = None, + prob: float = 1.0, ) -> None: self.id = id self.config = config self.network = network self.parent = parent + + if prob < 0 or 1 < prob: + raise ValueError("Probability must be a value in the range [0, 1]") + + self.prob = prob if children: self.children = children diff --git a/tests/functional/test_investment_pathway.py b/tests/functional/test_investment_pathway.py index 5774cae7..bbd25749 100644 --- a/tests/functional/test_investment_pathway.py +++ b/tests/functional/test_investment_pathway.py @@ -235,9 +235,9 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( This example models a case where investment decisions have to be made in 2030 and 2040. - In 2030, we have full knowledge of the existing assets - - In 2040, two equiprobable hypothesis are possible : - - A case where there is no change in the generation assets since 2030 (except te potential investment in 2030) - - A case where a base generation unit is present + - In 2040, two possible hypothesis are possible : + - P=0.2 => A case where there is no change in the generation assets since 2030 (except te potential investment in 2030) + - P=0.8 => A case where a base generation unit is present When taking the decision in 2030, we do not know which case will occur in 2040 and we seek the best decision given a risk criterion (the expectation here). @@ -292,15 +292,15 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( Case 1 : prob | investment | operational root: 1 x [ 100 x 100 + 10 x 300 ] - child A: + 0.5 x [ 50 x 100 + 10 x 400 (generator) + 5 x 200 (base)] - child B: + 0.5 x [ 50 x 100 + 10 x 400 (generator) + 10 000 x 200 (unsupplied energy)] - = 1 022 500 + child A: + 0.8 x [ 50 x 100 + 10 x 400 (generator) + 5 x 200 (base)] + child B: + 0.2 x [ 50 x 100 + 10 x 400 (generator) + 10 000 x 200 (unsupplied energy)] + = 422 800 Case 2 : prob | investment | operational root: 1 x [ 100 x 300 + 10 x 300 ] - child A: + 0.5 x [ 50 x 0 + 10 x 400 (generator) + 5 x 200 (base)] - child B: + 0.5 x [ 50 x 100 + 10 x 600 (generator)] - = 41 000 + child A: + 0.8 x [ 50 x 0 + 10 x 400 (generator) + 5 x 200 (base)] + child B: + 0.2 x [ 50 x 100 + 10 x 600 (generator)] + = 39 200 As investing less than 300 in the first stage would increase the unsupplied energy and lead to an increase in overall cost (-1 MW invested in 1st stage => + 1 MW unsupplied energy => +900/MW cost increase more or less), the optimal solution is to invest : @@ -443,12 +443,13 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( [TimeBlock(0, [0])], scenarios ) + # TODO Implement the prob behavior for the Expected Value dt_root = DecisionTreeNode("root", time_scenario_config, network_root) dt_child_A = DecisionTreeNode( - "childA", time_scenario_config, network_childA, parent=dt_root + "childA", time_scenario_config, network_childA, parent=dt_root, prob=0.8 ) dt_child_B = DecisionTreeNode( - "childB", time_scenario_config, network_childB, parent=dt_root + "childB", time_scenario_config, network_childB, parent=dt_root, prob=0.2 ) xpansion = build_benders_decomposed_problem( diff --git a/tests/unittests/test_generate_network_on_tree.py b/tests/unittests/test_generate_network_on_tree.py index 20d93a0b..c7d49f61 100644 --- a/tests/unittests/test_generate_network_on_tree.py +++ b/tests/unittests/test_generate_network_on_tree.py @@ -10,6 +10,7 @@ # # This file is part of the Antares project. +import pytest from andromede.simulation import TimeBlock from andromede.simulation.decision_tree import ( @@ -29,14 +30,26 @@ def test_generate_model_on_node() -> None: assert root.id == "root" assert root.parent is None + assert root.prob == 1.0 assert not root.children # No children - child = DecisionTreeNode("child", config, parent=root) + child = DecisionTreeNode("child", config, parent=root, prob=0.8) assert child.parent == root + assert child.prob == 0.8 assert child in root.children grandchild = DecisionTreeNode("grandchild", config, parent=child) assert grandchild.parent == child assert (grandchild not in root.children) and (grandchild in child.children) + + with pytest.raises(ValueError, match="Probability must be a value in the range"): + great_grandchild = DecisionTreeNode( + "greatgrandchild", config, parent=grandchild, prob=2.0 + ) + + with pytest.raises(ValueError, match="Probability must be a value in the range"): + great_grandchild = DecisionTreeNode( + "greatgrandchild", config, parent=grandchild, prob=-0.3 + ) From 7893e94a4570e0c945bc7443c717228683e22efe Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Thu, 20 Jun 2024 11:19:00 +0200 Subject: [PATCH 11/14] Remove test debug init --- tests/functional/test_investment_pathway.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/functional/test_investment_pathway.py b/tests/functional/test_investment_pathway.py index bbd25749..27133421 100644 --- a/tests/functional/test_investment_pathway.py +++ b/tests/functional/test_investment_pathway.py @@ -455,7 +455,6 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( xpansion = build_benders_decomposed_problem( dt_root, database, coupling_network=network_coupler ) - xpansion.initialise(is_debug=True) data = { "solution": { From 0f50a9bef0e0c56775f99c6c2d430f8e0fac967b Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Wed, 10 Jul 2024 15:19:14 +0200 Subject: [PATCH 12/14] Add prob to decision tree node --- src/andromede/simulation/decision_tree.py | 17 ++++++- .../test_generate_network_on_tree.py | 50 ++++++++++++++++++- 2 files changed, 64 insertions(+), 3 deletions(-) diff --git a/src/andromede/simulation/decision_tree.py b/src/andromede/simulation/decision_tree.py index f9404e5d..fc16e144 100644 --- a/src/andromede/simulation/decision_tree.py +++ b/src/andromede/simulation/decision_tree.py @@ -10,6 +10,7 @@ # # This file is part of the Antares project. +import math from dataclasses import dataclass from typing import Generator, Iterable, List, Optional @@ -48,7 +49,7 @@ def __init__( if prob < 0 or 1 < prob: raise ValueError("Probability must be a value in the range [0, 1]") - self.prob = prob + self.prob = prob * (parent.prob if parent is not None else 1) if children: self.children = children @@ -56,3 +57,17 @@ def traverse( self, depth: Optional[int] = None ) -> Generator["DecisionTreeNode", None, None]: yield from LevelOrderIter(self, maxlevel=depth) + + def is_leaves_prob_sum_one(self) -> bool: + if not self.children: + return True + + # Since we multiply the child's prob by the parent's prob + # in the constructor, the sum of the children prob should + # equal 1 * parent.prob if the values were set correctly + if not math.isclose(self.prob, sum(child.prob for child in self.children)): + return False + + # Recursively check if child nodes have their children's + # probability sum equal to one + return all(child.is_leaves_prob_sum_one() for child in self.children) diff --git a/tests/unittests/test_generate_network_on_tree.py b/tests/unittests/test_generate_network_on_tree.py index c7d49f61..6ac77b9b 100644 --- a/tests/unittests/test_generate_network_on_tree.py +++ b/tests/unittests/test_generate_network_on_tree.py @@ -20,7 +20,7 @@ from andromede.study.network import Network -def test_generate_model_on_node() -> None: +def test_decision_tree_generation() -> None: scenarios = 1 blocks = [TimeBlock(1, [0])] config = InterDecisionTimeScenarioConfig(blocks, scenarios) @@ -39,9 +39,10 @@ def test_generate_model_on_node() -> None: assert child.prob == 0.8 assert child in root.children - grandchild = DecisionTreeNode("grandchild", config, parent=child) + grandchild = DecisionTreeNode("grandchild", config, parent=child, prob=0.6) assert grandchild.parent == child + assert grandchild.prob == (0.8 * 0.6) assert (grandchild not in root.children) and (grandchild in child.children) with pytest.raises(ValueError, match="Probability must be a value in the range"): @@ -53,3 +54,48 @@ def test_generate_model_on_node() -> None: great_grandchild = DecisionTreeNode( "greatgrandchild", config, parent=grandchild, prob=-0.3 ) + + +def test_decision_tree_probabilities() -> None: + scenarios = 1 + blocks = [TimeBlock(1, [0])] + config = InterDecisionTimeScenarioConfig(blocks, scenarios) + network = Network("network_id") + + """ + root (p = 1) + |- l_child (p = 0.7) + | |- ll_child (p = 0.5) + | | `- lll_child (p = 1) + | | + | `- lr_child (p = 0.5) + | + `- r_child (p = 0.3) + |- rl_child (p = 0.4) + `- rr_child (p = 0.5) + """ + + # Root + root = DecisionTreeNode("root", config, network) + + # 1st level + l_child = DecisionTreeNode("l_child", config, parent=root, prob=0.7) + r_child = DecisionTreeNode("r_child", config, parent=root, prob=0.3) + + # 2nd level + ll_child = DecisionTreeNode("ll_child", config, parent=l_child, prob=0.5) + lr_child = DecisionTreeNode("lr_child", config, parent=l_child, prob=0.5) + + rl_child = DecisionTreeNode("rl_child", config, parent=r_child, prob=0.4) + rr_child = DecisionTreeNode("rr_child", config, parent=r_child, prob=0.5) + + # 3rd level + lll_child = DecisionTreeNode("lll_child", config, parent=ll_child, prob=1) + + assert ll_child.is_leaves_prob_sum_one() # One child with p = 1 + + assert l_child.is_leaves_prob_sum_one() # Two children w/ p1 = p2 = 0.5 + + assert not r_child.is_leaves_prob_sum_one() # Two children w/ p1 + p2 != 1 + + assert not root.is_leaves_prob_sum_one() From 9e92f5465f1cf10b57d3cd19ca91adf8626959fe Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Thu, 11 Jul 2024 13:21:57 +0200 Subject: [PATCH 13/14] Add ExpectedValueStrategy to Benders --- .../simulation/benders_decomposed.py | 7 ++++ src/andromede/simulation/optimization.py | 27 +++++++++++---- src/andromede/simulation/strategy.py | 33 ++++++++++++++++++- tests/functional/test_investment_pathway.py | 5 +-- 4 files changed, 61 insertions(+), 11 deletions(-) diff --git a/src/andromede/simulation/benders_decomposed.py b/src/andromede/simulation/benders_decomposed.py index 842dd10b..a7d6b1cd 100644 --- a/src/andromede/simulation/benders_decomposed.py +++ b/src/andromede/simulation/benders_decomposed.py @@ -32,6 +32,7 @@ ) from andromede.simulation.runner import BendersRunner, MergeMPSRunner from andromede.simulation.strategy import ( + ExpectedValue, InvestmentProblemStrategy, OperationalProblemStrategy, ) @@ -221,6 +222,9 @@ def build_benders_decomposed_problem( Returns a Benders Decomposed problem """ + if not decision_tree_root.is_leaves_prob_sum_one(): + raise RuntimeError("Decision tree must have leaves' probability sum equal one!") + null_time_block = TimeBlock( 0, [0] ) # Not necessary for master, but list must be non-empty @@ -234,6 +238,7 @@ def build_benders_decomposed_problem( problem_name="coupler", solver_id=solver_id, build_strategy=InvestmentProblemStrategy(), + risk_strategy=ExpectedValue(0.0), ) masters = [] # Benders Decomposed Master Problem @@ -252,6 +257,7 @@ def build_benders_decomposed_problem( solver_id=solver_id, build_strategy=InvestmentProblemStrategy(), decision_tree_node=tree_node.id, + risk_strategy=ExpectedValue(tree_node.prob), ) ) @@ -269,6 +275,7 @@ def build_benders_decomposed_problem( solver_id=solver_id, build_strategy=OperationalProblemStrategy(), decision_tree_node=tree_node.id, + risk_strategy=ExpectedValue(tree_node.prob), ) ) diff --git a/src/andromede/simulation/optimization.py b/src/andromede/simulation/optimization.py index bef4ad54..fc74376e 100644 --- a/src/andromede/simulation/optimization.py +++ b/src/andromede/simulation/optimization.py @@ -40,7 +40,12 @@ 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.strategy import ( + MergedProblemStrategy, + ModelSelectionStrategy, + RiskManagementStrategy, + UniformRisk, +) from andromede.simulation.time_block import TimeBlock from andromede.study.data import DataBase from andromede.study.network import Component, Network @@ -290,6 +295,7 @@ def __init__( scenarios: int, border_management: BlockBorderManagement, build_strategy: ModelSelectionStrategy = MergedProblemStrategy(), + risk_strategy: RiskManagementStrategy = UniformRisk(), decision_tree_node: str = "", ): self._network = network @@ -298,6 +304,7 @@ def __init__( self._scenarios = scenarios self._border_management = border_management self._build_strategy = build_strategy + self._risk_strategy = risk_strategy self._tree_node = decision_tree_node self._component_variables: Dict[TimestepComponentVariableKey, lp.Variable] = {} @@ -319,9 +326,13 @@ def tree_node(self) -> str: return self._tree_node @property - def strategy(self) -> ModelSelectionStrategy: + def build_strategy(self) -> ModelSelectionStrategy: return self._build_strategy + @property + def risk_strategy(self) -> RiskManagementStrategy: + return self._risk_strategy + def block_length(self) -> int: return len(self._block.timesteps) @@ -714,7 +725,7 @@ def _create_variables(self) -> None: component_context = self.context.create_component_context(component) model = component.model - for model_var in self.context.strategy.get_variables(model): + for model_var in self.context.build_strategy.get_variables(model): var_indexing = IndexingStructure( model_var.structure.time, model_var.structure.scenario ) @@ -770,7 +781,9 @@ def _create_variables(self) -> None: def _create_constraints(self) -> None: for component in self.context.network.all_components: - for constraint in self.context.strategy.get_constraints(component.model): + for constraint in self.context.build_strategy.get_constraints( + component.model + ): instantiated_expr = _instantiate_model_expression( constraint.expression, component.id, self.context ) @@ -798,14 +811,14 @@ def _create_objectives(self) -> None: component_context = self.context.create_component_context(component) model = component.model - for objective in self.context.strategy.get_objectives(model): + for objective in self.context.build_strategy.get_objectives(model): if objective is not None: _create_objective( self.solver, self.context, component, component_context, - objective, + self.context.risk_strategy(objective), ) def export_as_mps(self) -> str: @@ -825,6 +838,7 @@ def build_problem( border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, solver_id: str = "GLOP", build_strategy: ModelSelectionStrategy = MergedProblemStrategy(), + risk_strategy: RiskManagementStrategy = UniformRisk(), decision_tree_node: str = "", ) -> OptimizationProblem: """ @@ -841,6 +855,7 @@ def build_problem( scenarios, border_management, build_strategy, + risk_strategy, decision_tree_node, ) diff --git a/src/andromede/simulation/strategy.py b/src/andromede/simulation/strategy.py index 75e34c65..645ea646 100644 --- a/src/andromede/simulation/strategy.py +++ b/src/andromede/simulation/strategy.py @@ -13,7 +13,7 @@ from abc import ABC, abstractmethod from typing import Generator, Optional -from andromede.expression import ExpressionNode +from andromede.expression import ExpressionNode, literal from andromede.model import Constraint, Model, ProblemContext, Variable @@ -80,3 +80,34 @@ def get_objectives( self, model: Model ) -> Generator[Optional[ExpressionNode], None, None]: yield model.objective_operational_contribution + + +class RiskManagementStrategy(ABC): + """ + Abstract functor class for risk management + Its derived classes will implement risk measures: + - UniformRisk : The default case. All expressions have the same weight + - ExpectedValue : Computes the product prob * expression + TODO For now, it will only take into account the Expected Value + TODO In the future could have other risk measures? + """ + + def __call__(self, expr: ExpressionNode) -> ExpressionNode: + return self._modify_expression(expr) + + @abstractmethod + def _modify_expression(self, expr: ExpressionNode) -> ExpressionNode: + ... + + +class UniformRisk(RiskManagementStrategy): + def _modify_expression(self, expr: ExpressionNode) -> ExpressionNode: + return expr + + +class ExpectedValue(RiskManagementStrategy): + def __init__(self, prob: float) -> None: + self._prob = prob + + def _modify_expression(self, expr: ExpressionNode) -> ExpressionNode: + return literal(self._prob) * expr diff --git a/tests/functional/test_investment_pathway.py b/tests/functional/test_investment_pathway.py index 27133421..dd77f5e2 100644 --- a/tests/functional/test_investment_pathway.py +++ b/tests/functional/test_investment_pathway.py @@ -309,9 +309,6 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( - 100 in child B """ - # Either we duplicate all network for each node : Lots of duplications - # or we index all data, parameters, variables by the resolution node : Make the data structure dependent of the resolution tree... - database = DataBase() database.add_data("N", "spillage_cost", ConstantData(10)) database.add_data("N", "ens_cost", ConstantData(10_000)) @@ -458,7 +455,7 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( data = { "solution": { - "overall_cost": 49_000, + "overall_cost": 39_200, "values": { "root_CAND_delta_invest": 300, "childA_CAND_delta_invest": 0, From 244fa11d676b357c1c5c580e1208810de7e7079b Mon Sep 17 00:00:00 2001 From: Ian Menezes Date: Tue, 16 Jul 2024 11:20:18 +0200 Subject: [PATCH 14/14] Retour of PR comments --- src/andromede/expression/expression.py | 4 ++ src/andromede/model/constraint.py | 7 +++ .../simulation/benders_decomposed.py | 14 ++++-- src/andromede/simulation/optimization.py | 18 ++++---- tests/functional/test_andromede.py | 2 - tests/functional/test_investment_pathway.py | 44 +++++++++++++++++-- 6 files changed, 72 insertions(+), 17 deletions(-) diff --git a/src/andromede/expression/expression.py b/src/andromede/expression/expression.py index 4aa013f1..7262326d 100644 --- a/src/andromede/expression/expression.py +++ b/src/andromede/expression/expression.py @@ -226,6 +226,10 @@ def is_unbound(expr: ExpressionNode) -> bool: return isinstance(expr, LiteralNode) and (abs(expr.value) == float("inf")) +def is_non_negative(expr: ExpressionNode) -> bool: + return isinstance(expr, LiteralNode) and (expr.value >= 0) + + @dataclass(frozen=True, eq=False) class UnaryOperatorNode(ExpressionNode): operand: ExpressionNode diff --git a/src/andromede/model/constraint.py b/src/andromede/model/constraint.py index ea87116c..0f33aa7f 100644 --- a/src/andromede/model/constraint.py +++ b/src/andromede/model/constraint.py @@ -18,6 +18,7 @@ Comparator, ComparisonNode, ExpressionNode, + is_non_negative, is_unbound, literal, ) @@ -69,5 +70,11 @@ def __post_init__( f"The bounds of a constraint should not contain variables, {print_expr(bound)} was given." ) + if is_unbound(self.lower_bound) and is_non_negative(self.lower_bound): + raise ValueError("Lower bound should not be +Inf") + + if is_unbound(self.upper_bound) and not is_non_negative(self.upper_bound): + raise ValueError("Upper bound should not be -Inf") + def replicate(self, /, **changes: Any) -> "Constraint": return replace(self, **changes) diff --git a/src/andromede/simulation/benders_decomposed.py b/src/andromede/simulation/benders_decomposed.py index a7d6b1cd..57a03ff4 100644 --- a/src/andromede/simulation/benders_decomposed.py +++ b/src/andromede/simulation/benders_decomposed.py @@ -52,6 +52,7 @@ class BendersDecomposedProblem: emplacement: pathlib.Path output_path: pathlib.Path + structure_filename: str solution: Optional[BendersSolution] is_merged: bool @@ -62,12 +63,14 @@ def __init__( subproblems: List[OptimizationProblem], emplacement: str = "outputs/lp", output_path: str = "expansion", + struct_filename: str = "structure.txt", ) -> None: self.master = master self.subproblems = subproblems self.emplacement = pathlib.Path(emplacement) self.output_path = pathlib.Path(output_path) + self.structure_filename = struct_filename self.solution = None self.is_merged = False @@ -124,7 +127,7 @@ def export_options( "SLAVE_WEIGHT": "CONSTANT", "SLAVE_WEIGHT_VALUE": 1, "MASTER_NAME": f"{self.master.name}", - "STRUCTURE_FILE": "structure.txt", + "STRUCTURE_FILE": f"{self.structure_filename}", "INPUTROOT": ".", "CSV_NAME": "benders_output_trace", "BOUND_ALPHA": True, @@ -155,7 +158,9 @@ def initialise( serialize( f"{subproblem.name}.mps", subproblem.export_as_mps(), self.emplacement ) - serialize(f"structure.txt", self.export_structure(), self.emplacement) + serialize( + f"{self.structure_filename}", self.export_structure(), self.emplacement + ) serialize_json( "options.json", self.export_options(solver_name=solver_name, log_level=log_level), @@ -215,6 +220,7 @@ def build_benders_decomposed_problem( border_management: BlockBorderManagement = BlockBorderManagement.CYCLE, solver_id: str = "GLOP", coupling_network: Network = Network(""), + struct_filename: str = "structure.txt", ) -> BendersDecomposedProblem: """ Entry point to build the xpansion problem for a time period @@ -281,4 +287,6 @@ def build_benders_decomposed_problem( master = fusion_problems(masters, coupler) - return BendersDecomposedProblem(master, subproblems) + return BendersDecomposedProblem( + master, subproblems, struct_filename=struct_filename + ) diff --git a/src/andromede/simulation/optimization.py b/src/andromede/simulation/optimization.py index fc74376e..b53d71c7 100644 --- a/src/andromede/simulation/optimization.py +++ b/src/andromede/simulation/optimization.py @@ -742,17 +742,17 @@ def _create_variables(self) -> None: # Set solver var name # Externally, for the Solver, this variable will have a full name - # Internally, it will be indexed by a structure that into account + # Internally, it will be indexed by a structure that takes into account # the component id, variable name, timestep and scenario separately solver_var_name: str = f"{model_var.name}" if component.id: - solver_var_name = f"{component.id}_" + solver_var_name + solver_var_name = f"{component.id}_{solver_var_name}" if self.context.tree_node: - solver_var_name = f"{self.context.tree_node}_" + solver_var_name + solver_var_name = f"{self.context.tree_node}_{solver_var_name}" for block_timestep in self.context.get_time_indices(var_indexing): if self.context.block_length() > 1: - solver_var_name = solver_var_name + f"_t{block_timestep}" + solver_var_name = f"{solver_var_name}_t{block_timestep}" for scenario in self.context.get_scenario_indices(var_indexing): lower_bound = -self.solver.infinity() @@ -767,7 +767,7 @@ def _create_variables(self) -> None: ).get_value(block_timestep, scenario) if self.context.scenarios > 1: - solver_var_name = solver_var_name + f"_s{scenario}" + solver_var_name = f"{solver_var_name}_s{scenario}" # TODO: Add BoolVar or IntVar if the variable is specified to be integer or bool solver_var = self.solver.NumVar( @@ -873,7 +873,7 @@ def fusion_problems( root_master.name = "master" root_vars: Dict[str, lp.Variable] = dict() - root_cstrs: Dict[str, lp.Constraint] = dict() + root_constraints: Dict[str, lp.Constraint] = dict() root_objective = root_master.solver.Objective() # We stock the coupler's variables to check for @@ -910,11 +910,11 @@ def fusion_problems( # If variable present in constraint, we add the constraint to root if coeff != 0: key = f"{master.name}_{cstr.name()}" - if key not in root_cstrs: - root_cstrs[key] = root_master.solver.Constraint( + if key not in root_constraints: + root_constraints[key] = root_master.solver.Constraint( cstr.Lb(), cstr.Ub(), key ) - root_cstr = root_cstrs[key] + root_cstr = root_constraints[key] root_cstr.SetCoefficient(root_var, coeff) obj_coeff = objective.GetCoefficient(var) diff --git a/tests/functional/test_andromede.py b/tests/functional/test_andromede.py index ba033a32..300ff88d 100644 --- a/tests/functional/test_andromede.py +++ b/tests/functional/test_andromede.py @@ -30,8 +30,6 @@ from andromede.model.model import PortFieldDefinition, PortFieldId from andromede.simulation import ( BlockBorderManagement, - DecisionTreeNode, - InterDecisionTimeScenarioConfig, OutputValues, TimeBlock, build_problem, diff --git a/tests/functional/test_investment_pathway.py b/tests/functional/test_investment_pathway.py index dd77f5e2..d2a25953 100644 --- a/tests/functional/test_investment_pathway.py +++ b/tests/functional/test_investment_pathway.py @@ -11,7 +11,6 @@ # This file is part of the Antares project. import pytest -from anytree import Node as TreeNode from andromede.expression import literal, var from andromede.expression.indexing_structure import IndexingStructure @@ -72,6 +71,40 @@ def test_investment_pathway_on_sequential_nodes( demand: Component, candidate: Component, ) -> None: + """ + A first simple test on the investment pathway + Here, only two nodes are represented, a parent and a child nodes + with probability one of going from parent to child + + The goal here is to show that, in the parent node, the demand is already met + by the existing fixed production. However, for the child node without any new + investment, it would create some unsupplied energy, which is very expensive. + + The investment on the child node, even though enough for the demand, is also + more expensive than on the parent (this should represent a late investment fee). + + To minimize the expected cost in this 2-node tree, one should expect the maximum + investment on the parent node, and the rest on the child node. + + Here below the values used: + + PARENT | CHILD + Demand (MW): 100 200 + Fixed prod (MW): 100 100 + Max invest (MW): 80 100 + Op cost ($): 10 10 + Invest cost ($): 100 300 + ENS cost ($): 10000 10000 + + The solution should be: + + prob | investment | operational + parent: 1 x [ 100 x 80 + 10 x 100 ] + child : + 1 x [ 300 x 20 + 10 x 200 ] + = 17 000 + + """ + # === Populating Database === database = DataBase() database.add_data("N", "spillage_cost", ConstantData(10)) database.add_data("N", "ens_cost", ConstantData(10000)) @@ -112,6 +145,8 @@ def test_investment_pathway_on_sequential_nodes( ), ) + # === Coupling model === + # Used between nodes in the decision tree COUPLING_MODEL = model( id="COUPLING", variables=[ @@ -156,6 +191,7 @@ def test_investment_pathway_on_sequential_nodes( ], ) + # === Network === network_coupling = Network("coupling_test") network_coupling.add_component(create_component(model=COUPLING_MODEL, id="")) @@ -187,6 +223,7 @@ def test_investment_pathway_on_sequential_nodes( PortRef(candidate_chd, "balance_port"), PortRef(node, "balance_port") ) + # === Decision tree creation === config = InterDecisionTimeScenarioConfig([TimeBlock(0, [0])], 1) decision_tree_par = DecisionTreeNode("parent", config, network_par) @@ -194,6 +231,7 @@ def test_investment_pathway_on_sequential_nodes( "child", config, network_chd, parent=decision_tree_par ) + # === Build problem === xpansion = build_benders_decomposed_problem( decision_tree_par, database, coupling_network=network_coupling ) @@ -211,6 +249,7 @@ def test_investment_pathway_on_sequential_nodes( } solution = BendersSolution(data) + # === Run === assert xpansion.run() decomposed_solution = xpansion.solution if decomposed_solution is not None: # For mypy only @@ -236,7 +275,7 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( This example models a case where investment decisions have to be made in 2030 and 2040. - In 2030, we have full knowledge of the existing assets - In 2040, two possible hypothesis are possible : - - P=0.2 => A case where there is no change in the generation assets since 2030 (except te potential investment in 2030) + - P=0.2 => A case where there is no change in the generation assets since 2030 (except the potential investment in 2030) - P=0.8 => A case where a base generation unit is present When taking the decision in 2030, we do not know which case will occur in 2040 @@ -440,7 +479,6 @@ def test_investment_pathway_on_a_tree_with_one_root_two_children( [TimeBlock(0, [0])], scenarios ) - # TODO Implement the prob behavior for the Expected Value dt_root = DecisionTreeNode("root", time_scenario_config, network_root) dt_child_A = DecisionTreeNode( "childA", time_scenario_config, network_childA, parent=dt_root, prob=0.8