Benders decomposition for investment problems (#12)

- Add benders-decomposed problem resolution, which will call
  benders tool from antares-xpansion.
- Variables, constraints and objectives can be tagged as investment
  or operational.
- Problem building from the data can now be customized with a strategy to
  select which variables etc. to include. This strategy has 2
  implementations: operational and investment.

.gitignore

@@ -8,3 +8,5 @@ venv
 .env
 .coverage
 coverage.xml
+bin
+outputs

src/andromede/libs/standard.py

@@ -41,6 +41,28 @@
     ],
 )
 
+NODE_WITH_SPILL_AND_ENS_MODEL = model(
+    id="NODE_WITH_SPILL_AND_ENS_MODEL",
+    parameters=[float_parameter("spillage_cost"), float_parameter("ens_cost")],
+    variables=[
+        float_variable("spillage", lower_bound=literal(0)),
+        float_variable("unsupplied_energy", lower_bound=literal(0)),
+    ],
+    ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")],
+    binding_constraints=[
+        Constraint(
+            name="Balance",
+            expression=port_field("balance_port", "flow").sum_connections()
+            == var("spillage") - var("unsupplied_energy"),
+        )
+    ],
+    objective_operational_contribution=(
+        param("spillage_cost") * var("spillage")
+        + param("ens_cost") * var("unsupplied_energy")
+    )
+    .sum()
+    .expec(),
+)
 """
 A standard model for a linear cost generation, limited by a maximum generation.
 """
@@ -63,7 +85,9 @@
             name="Max generation", expression=var("generation") <= param("p_max")
         ),
     ],
-    objective_contribution=(param("cost") * var("generation")).sum().expec(),
+    objective_operational_contribution=(param("cost") * var("generation"))
+    .sum()
+    .expec(),
 )
 
 """
@@ -133,7 +157,9 @@
             lower_bound=literal(0),
         ),  # To test both ways of setting constraints
     ],
-    objective_contribution=(param("cost") * var("generation")).sum().expec(),
+    objective_operational_contribution=(param("cost") * var("generation"))
+    .sum()
+    .expec(),
 )
 
 # For now, no starting cost
@@ -159,17 +185,17 @@
             "nb_on",
             lower_bound=literal(0),
             upper_bound=param("nb_units_max"),
-            structural_type=ANTICIPATIVE_TIME_VARYING,
+            structure=ANTICIPATIVE_TIME_VARYING,
         ),
         int_variable(
             "nb_stop",
             lower_bound=literal(0),
-            structural_type=ANTICIPATIVE_TIME_VARYING,
+            structure=ANTICIPATIVE_TIME_VARYING,
         ),
         int_variable(
             "nb_start",
             lower_bound=literal(0),
-            structural_type=ANTICIPATIVE_TIME_VARYING,
+            structure=ANTICIPATIVE_TIME_VARYING,
         ),
     ],
     ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")],
@@ -208,7 +234,9 @@
         )
         # It also works by writing ExpressionRange(-param("d_min_down") + 1, 0) as ExpressionRange's __post_init__ wraps integers to literal nodes. However, MyPy does not seem to infer that ExpressionRange's attributes are necessarily of ExpressionNode type and raises an error if the arguments in the constructor are integer (whereas it runs correctly), this why we specify it here with literal(0) instead of 0.
     ],
-    objective_contribution=(param("cost") * var("generation")).sum().expec(),
+    objective_operational_contribution=(param("cost") * var("generation"))
+    .sum()
+    .expec(),
 )
 
 # Same model as previous one, except that starting/stopping variables are now non anticipative
@@ -234,17 +262,17 @@
             "nb_on",
             lower_bound=literal(0),
             upper_bound=param("nb_units_max"),
-            structural_type=NON_ANTICIPATIVE_TIME_VARYING,
+            structure=NON_ANTICIPATIVE_TIME_VARYING,
         ),
         int_variable(
             "nb_stop",
             lower_bound=literal(0),
-            structural_type=NON_ANTICIPATIVE_TIME_VARYING,
+            structure=NON_ANTICIPATIVE_TIME_VARYING,
         ),
         int_variable(
             "nb_start",
             lower_bound=literal(0),
-            structural_type=NON_ANTICIPATIVE_TIME_VARYING,
+            structure=NON_ANTICIPATIVE_TIME_VARYING,
         ),
     ],
     ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")],
@@ -282,7 +310,9 @@
             <= param("nb_units_max").shift(-param("d_min_down")) - var("nb_on"),
         ),
     ],
-    objective_contribution=(param("cost") * var("generation")).sum().expec(),
+    objective_operational_contribution=(param("cost") * var("generation"))
+    .sum()
+    .expec(),
 )
 
 SPILLAGE_MODEL = model(
@@ -296,7 +326,7 @@
             definition=-var("spillage"),
         )
     ],
-    objective_contribution=(param("cost") * var("spillage")).sum().expec(),
+    objective_operational_contribution=(param("cost") * var("spillage")).sum().expec(),
 )
 
 UNSUPPLIED_ENERGY_MODEL = model(
@@ -310,7 +340,9 @@
             definition=var("unsupplied_energy"),
         )
     ],
-    objective_contribution=(param("cost") * var("unsupplied_energy")).sum().expec(),
+    objective_operational_contribution=(param("cost") * var("unsupplied_energy"))
+    .sum()
+    .expec(),
 )
 
 # Simplified model
@@ -357,5 +389,5 @@
             == param("inflows"),
         ),
     ],
-    objective_contribution=literal(0),  # Implcitement nul ?
+    objective_operational_contribution=literal(0),  # Implcitement nul ?
 )

src/andromede/model/__init__.py

@@ -10,8 +10,9 @@
 #
 # This file is part of the Antares project.
 
+from .common import ProblemContext, ValueType
 from .constraint import Constraint
 from .model import Model, ModelPort, model
-from .parameter import Parameter, ParameterValueType, float_parameter, int_parameter
+from .parameter import Parameter, float_parameter, int_parameter
 from .port import PortField, PortType
-from .variable import Variable, VariableValueType, float_variable, int_variable
+from .variable import Variable, float_variable, int_variable

src/andromede/model/common.py

@@ -0,0 +1,28 @@
+# Copyright (c) 2024, RTE (https://www.rte-france.com)
+#
+# See AUTHORS.txt
+#
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#
+# SPDX-License-Identifier: MPL-2.0
+#
+# This file is part of the Antares project.
+
+"""
+Module for common classes used in models.
+"""
+from enum import Enum
+
+
+class ValueType(Enum):
+    FLOAT = "FLOAT"
+    INTEGER = "INTEGER"
+    # Needs more ?
+
+
+class ProblemContext(Enum):
+    OPERATIONAL = 0
+    INVESTMENT = 1
+    COUPLING = 2

src/andromede/model/constraint.py

@@ -20,6 +20,7 @@
     literal,
 )
 from andromede.expression.print import print_expr
+from andromede.model.common import ProblemContext
 
 
 class Constraint:
@@ -33,15 +34,19 @@ class Constraint:
     expression: ExpressionNode
     lower_bound: ExpressionNode
     upper_bound: ExpressionNode
+    context: ProblemContext
 
     def __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:
                 raise ValueError(

src/andromede/model/model.py

@@ -78,6 +78,23 @@ def get_component_variable_structure(
     return Provider()
 
 
+def _is_objective_contribution_valid(
+    model: "Model", objective_contribution: ExpressionNode
+) -> bool:
+    if not is_linear(objective_contribution):
+        raise ValueError("Objective contribution must be a linear expression.")
+
+    data_structure_provider = _make_structure_provider(model)
+    objective_structure = compute_indexation(
+        objective_contribution, data_structure_provider
+    )
+
+    if objective_structure != IndexingStructure(time=False, scenario=False):
+        raise ValueError("Objective contribution should be a real-valued expression.")
+    # TODO: We should also check that the number of instances is equal to 1, but this would require a linearization here, do not want to do that for now...
+    return True
+
+
 @dataclass(frozen=True)
 class ModelPort:
     """
@@ -129,26 +146,23 @@ class Model:
     inter_block_dyn: bool = False
     parameters: Dict[str, Parameter] = field(default_factory=dict)
     variables: Dict[str, Variable] = field(default_factory=dict)
-    objective_contribution: Optional[ExpressionNode] = None
+    objective_operational_contribution: Optional[ExpressionNode] = None
+    objective_investment_contribution: Optional[ExpressionNode] = None
     ports: Dict[str, ModelPort] = field(default_factory=dict)  # key = port name
     port_fields_definitions: Dict[PortFieldId, PortFieldDefinition] = field(
         default_factory=dict
     )
 
     def __post_init__(self) -> None:
-        if self.objective_contribution:
-            if not is_linear(self.objective_contribution):
-                raise ValueError("Objective contribution must be a linear expression.")
+        if self.objective_operational_contribution:
+            _is_objective_contribution_valid(
+                self, self.objective_operational_contribution
+            )
 
-            data_structure_provider = _make_structure_provider(self)
-            objective_structure = compute_indexation(
-                self.objective_contribution, data_structure_provider
+        if self.objective_investment_contribution:
+            _is_objective_contribution_valid(
+                self, self.objective_investment_contribution
             )
-            if objective_structure != IndexingStructure(time=False, scenario=False):
-                raise ValueError(
-                    "Objective contribution should be a real-valued expression."
-                )
-            # TODO: We should also check that the number of instances is equal to 1, but this would require a linearization here, do not want to do that for now...
 
         for definition in self.port_fields_definitions.values():
             port_name = definition.port_field.port_name
@@ -176,7 +190,8 @@ def model(
     binding_constraints: Optional[Iterable[Constraint]] = None,
     parameters: Optional[Iterable[Parameter]] = None,
     variables: Optional[Iterable[Variable]] = None,
-    objective_contribution: Optional[ExpressionNode] = None,
+    objective_operational_contribution: Optional[ExpressionNode] = None,
+    objective_investment_contribution: Optional[ExpressionNode] = None,
     inter_block_dyn: bool = False,
     ports: Optional[Iterable[ModelPort]] = None,
     port_fields_definitions: Optional[Iterable[PortFieldDefinition]] = None,
@@ -202,7 +217,8 @@ def model(
         else {},
         parameters={p.name: p for p in parameters} if parameters else {},
         variables={v.name: v for v in variables} if variables else {},
-        objective_contribution=objective_contribution,
+        objective_operational_contribution=objective_operational_contribution,
+        objective_investment_contribution=objective_investment_contribution,
         inter_block_dyn=inter_block_dyn,
         ports=existing_port_names,
         port_fields_definitions={d.port_field: d for d in port_fields_definitions}

src/andromede/model/parameter.py

@@ -11,15 +11,9 @@
 # This file is part of the Antares project.
 
 from dataclasses import dataclass
-from enum import Enum
 
 from andromede.expression.indexing_structure import IndexingStructure
-
-
-class ParameterValueType(Enum):
-    FLOAT = "FLOAT"
-    INTEGER = "INTEGER"
-    # Needs more ?
+from andromede.model.common import ValueType
 
 
 @dataclass(frozen=True)
@@ -31,19 +25,19 @@ class Parameter:
     """
 
     name: str
-    type: ParameterValueType
+    type: ValueType
     structure: IndexingStructure
 
 
 def int_parameter(
     name: str,
     structure: IndexingStructure = IndexingStructure(True, True),
 ) -> Parameter:
-    return Parameter(name, ParameterValueType.INTEGER, structure)
+    return Parameter(name, ValueType.INTEGER, structure)
 
 
 def float_parameter(
     name: str,
     structure: IndexingStructure = IndexingStructure(True, True),
 ) -> Parameter:
-    return Parameter(name, ParameterValueType.FLOAT, structure)
+    return Parameter(name, ValueType.FLOAT, structure)