Merge branch 'main' into feature/executable

src/andromede/expression/indexing_structure.py

@@ -26,3 +26,15 @@ def __or__(self, other: "IndexingStructure") -> "IndexingStructure":
         time = self.time or other.time
         scenario = self.scenario or other.scenario
         return IndexingStructure(time, scenario)
+
+    def is_time_varying(self) -> bool:
+        return self.time
+
+    def is_scenario_varying(self) -> bool:
+        return self.scenario
+
+    def is_time_scenario_varying(self) -> bool:
+        return self.is_time_varying() and self.is_scenario_varying()
+
+    def is_constant(self) -> bool:
+        return (not self.is_time_varying()) and (not self.is_scenario_varying())

src/andromede/libs/standard_sc.py

@@ -211,3 +211,81 @@
         )
     ],
 )
+
+SHORT_TERM_STORAGE_COMPLEX = model(
+    id="STS_COMPLEX",
+    parameters=[
+        float_parameter("p_max_injection"),
+        float_parameter("p_max_withdrawal"),
+        float_parameter("level_min"),
+        float_parameter("level_max"),
+        float_parameter("inflows"),
+        float_parameter(
+            "efficiency"
+        ),  # Should be constant, but time-dependent values should work as well
+        float_parameter("withdrawal_penality"),
+        float_parameter("level_penality"),
+        float_parameter("Pgrad+i_penality"),
+        float_parameter("Pgrad-i_penality"),
+        float_parameter("Pgrad+s_penality"),
+        float_parameter("Pgrad-s_penality"),
+    ],
+    variables=[
+        float_variable(
+            "injection", lower_bound=literal(0), upper_bound=param("p_max_injection")
+        ),
+        float_variable(
+            "withdrawal", lower_bound=literal(0), upper_bound=param("p_max_withdrawal")
+        ),
+        float_variable(
+            "level", lower_bound=param("level_min"), upper_bound=param("level_max")
+        ),
+        float_variable("Pgrad+i", lower_bound=literal(0)),
+        float_variable("Pgrad-i", lower_bound=literal(0)),
+        float_variable("Pgrad+s", lower_bound=literal(0)),
+        float_variable("Pgrad-s", lower_bound=literal(0)),
+    ],
+    ports=[ModelPort(port_type=BALANCE_PORT_TYPE, port_name="balance_port")],
+    port_fields_definitions=[
+        PortFieldDefinition(
+            port_field=PortFieldId("balance_port", "flow"),
+            definition=var("withdrawal") - var("injection"),
+        )
+    ],
+    constraints=[
+        Constraint(
+            name="Level",
+            expression=var("level")
+            - var("level").shift(-1)
+            - param("efficiency") * var("injection")
+            + var("withdrawal")
+            == param("inflows"),
+        ),
+        Constraint(
+            "Pgrad+i min",
+            var("Pgrad+i") >= var("injection") - var("injection").shift(-1),
+        ),
+        Constraint(
+            "Pgrad-i min",
+            var("Pgrad-i") >= var("injection").shift(-1) - var("injection"),
+        ),
+        Constraint(
+            "Pgrad+s min",
+            var("Pgrad+s") >= var("withdrawal") - var("withdrawal").shift(-1),
+        ),
+        Constraint(
+            "Pgrad-s min",
+            var("Pgrad-s") >= var("withdrawal").shift(-1) - var("withdrawal"),
+        ),
+    ],
+    objective_operational_contribution=(
+        param("level_penality") * var("level")
+        + param("withdrawal_penality") * var("withdrawal")
+        + param("Pgrad+i_penality") * var("Pgrad+i")
+        + param("Pgrad-i_penality") * var("Pgrad-i")
+        + param("Pgrad+s_penality") * var("Pgrad+s")
+        + param("Pgrad-s_penality") * var("Pgrad-s")
+    )
+    .sum()
+    .expec(),
+)

src/andromede/model/common.py

@@ -19,7 +19,7 @@
 class ValueType(Enum):
     FLOAT = "FLOAT"
     INTEGER = "INTEGER"
-    # Needs more ?
+    BOOL = "BOOL"
 
 
 class ProblemContext(Enum):

src/andromede/model/variable.py

@@ -13,12 +13,9 @@
 from dataclasses import dataclass
 from typing import Any, Optional
 
-from andromede.expression import ExpressionNode
+from andromede.expression import ExpressionNode, literal
 from andromede.expression.degree import is_constant
-from andromede.expression.equality import (
-    expressions_equal,
-    expressions_equal_if_present,
-)
+from andromede.expression.equality import expressions_equal_if_present
 from andromede.expression.indexing_structure import IndexingStructure
 from andromede.model.common import ProblemContext, ValueType
 
@@ -66,6 +63,14 @@ def int_variable(
     )
 
 
+def bool_var(
+    name: str,
+    structure: IndexingStructure = IndexingStructure(True, True),
+    context: ProblemContext = ProblemContext.OPERATIONAL,
+) -> Variable:
+    return Variable(name, ValueType.BOOL, literal(0), literal(1), structure, context)
+
+
 def float_variable(
     name: str,
     lower_bound: Optional[ExpressionNode] = None,

src/andromede/simulation/optimization.py

@@ -15,10 +15,11 @@
 into a mathematical optimization problem.
 """
 
+import math
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
 from enum import Enum
-from typing import Dict, Iterable, List, Optional, Type
+from typing import Dict, Iterable, List, Optional
 
 import ortools.linear_solver.pywraplp as lp
 
@@ -721,11 +722,10 @@ def _create_variables(self) -> None:
             model = component.model
 
             for model_var in self.strategy.get_variables(model):
-                var_indexing = IndexingStructure(
-                    model_var.structure.time, model_var.structure.scenario
-                )
+                var_indexing = model_var.structure
                 instantiated_lb_expr = None
                 instantiated_ub_expr = None
+
                 if model_var.lower_bound:
                     instantiated_lb_expr = _instantiate_model_expression(
                         model_var.lower_bound, component.id, self.context
@@ -734,7 +734,18 @@ def _create_variables(self) -> None:
                     instantiated_ub_expr = _instantiate_model_expression(
                         model_var.upper_bound, component.id, self.context
                     )
+
+                var_name: str = f"{model_var.name}"
+                component_prefix = f"{component.id}_" if component.id else ""
+
                 for block_timestep in self.context.get_time_indices(var_indexing):
+                    block_suffix = (
+                        f"_t{block_timestep}"
+                        if var_indexing.is_time_varying()
+                        and (self.context.block_length() > 1)
+                        else ""
+                    )
+
                     for scenario in self.context.get_scenario_indices(var_indexing):
                         lower_bound = -self.solver.infinity()
                         upper_bound = self.solver.infinity()
@@ -747,23 +758,46 @@ def _create_variables(self) -> None:
                                 instantiated_ub_expr
                             ).get_value(block_timestep, scenario)
 
-                        # TODO: Add BoolVar or IntVar if the variable is specified to be integer or bool
+                        scenario_suffix = (
+                            f"_s{scenario}"
+                            if var_indexing.is_scenario_varying()
+                            and (self.context.scenarios > 1)
+                            else ""
+                        )
+
+                        # 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 = None
-                        if model_var.data_type == ValueType.FLOAT:
-                            solver_var = self.solver.NumVar(
-                                lower_bound,
-                                upper_bound,
-                                f"{component.id}_{model_var.name}_t{block_timestep}_s{scenario}",
+                        solver_var_name = f"{component_prefix}{var_name}{block_suffix}{scenario_suffix}"
+
+                        if math.isclose(lower_bound, upper_bound):
+                            raise ValueError(
+                                f"Upper and lower bounds of variable {solver_var_name} have the same value: {lower_bound}"
+                            )
+                        elif lower_bound > upper_bound:
+                            raise ValueError(
+                                f"Upper bound ({upper_bound}) must be strictly greater than lower bound ({lower_bound}) for variable {solver_var_name}"
+                            )
+
+                        if model_var.data_type == ValueType.BOOL:
+                            solver_var = self.solver.BoolVar(
+                                solver_var_name,
                             )
                         elif model_var.data_type == ValueType.INTEGER:
                             solver_var = self.solver.IntVar(
                                 lower_bound,
                                 upper_bound,
-                                f"{component.id}_{model_var.name}_t{block_timestep}_s{scenario}",
+                                solver_var_name,
                             )
+                        else:
+                            solver_var = self.solver.NumVar(
+                                lower_bound,
+                                upper_bound,
+                                solver_var_name,
+                            )
+
                         component_context.add_variable(
                             block_timestep, scenario, model_var.name, solver_var
                         )

tests/functional/conftest.py

@@ -14,7 +14,7 @@
 import pytest
 
 from andromede.model.parsing import parse_yaml_library
-from andromede.model.resolve_library import resolve_library
+from andromede.model.resolve_library import Library, resolve_library
 
 
 @pytest.fixture(scope="session")
@@ -23,7 +23,7 @@ def libs_dir() -> Path:
 
 
 @pytest.fixture(scope="session")
-def lib(libs_dir: Path):
+def lib(libs_dir: Path) -> Library:
     lib_file = libs_dir / "lib.yml"
 
     with lib_file.open() as f:

tests/functional/test_andromede.py

@@ -173,6 +173,22 @@ def test_variable_bound() -> None:
     status = problem.solver.Solve()
     assert status == problem.solver.INFEASIBLE  # Infeasible
 
+    network = create_one_node_network(generator_model)
+    database = create_simple_database(max_generation=0)  # Equal upper and lower bounds
+    with pytest.raises(
+        ValueError,
+        match="Upper and lower bounds of variable G_generation have the same value: 0",
+    ):
+        problem = build_problem(network, database, TimeBlock(1, [0]), 1)
+
+    network = create_one_node_network(generator_model)
+    database = create_simple_database(max_generation=-10)
+    with pytest.raises(
+        ValueError,
+        match=r"Upper bound \(-10\) must be strictly greater than lower bound \(0\) for variable G_generation",
+    ):
+        problem = build_problem(network, database, TimeBlock(1, [0]), 1)
+
 
 def generate_data(
     efficiency: float, horizon: int, scenarios: int
@@ -181,7 +197,7 @@ def generate_data(
     for scenario in range(scenarios):
         for absolute_timestep in range(horizon):
             if absolute_timestep == 0:
-                data[TimeScenarioIndex(absolute_timestep, scenario)] = -18
+                data[TimeScenarioIndex(absolute_timestep, scenario)] = -18.0
             else:
                 data[TimeScenarioIndex(absolute_timestep, scenario)] = 2 * efficiency
 

tests/functional/test_andromede_yml.py

@@ -4,6 +4,7 @@
 from andromede.expression import literal, param, var
 from andromede.expression.indexing_structure import IndexingStructure
 from andromede.model import Model, ModelPort, float_parameter, float_variable, model
+from andromede.model.library import Library
 from andromede.model.model import PortFieldDefinition, PortFieldId
 from andromede.simulation import (
     BlockBorderManagement,
@@ -23,7 +24,7 @@
 )
 
 
-def test_network(lib) -> None:
+def test_network(lib: Library) -> None:
     network = Network("test")
     assert network.id == "test"
     assert list(network.nodes) == []
@@ -47,7 +48,7 @@ def test_network(lib) -> None:
         network.get_component("unknown")
 
 
-def test_basic_balance(lib) -> None:
+def test_basic_balance(lib: Library) -> None:
     """
     Balance on one node with one fixed demand and one generation, on 1 timestep.
     """
@@ -88,7 +89,7 @@ def test_basic_balance(lib) -> None:
     assert problem.solver.Objective().Value() == 3000
 
 
-def test_link(lib) -> None:
+def test_link(lib: Library) -> None:
     """
     Balance on one node with one fixed demand and one generation, on 1 timestep.
     """
@@ -146,7 +147,7 @@ def test_link(lib) -> None:
             assert variable.solution_value() == -100
 
 
-def test_stacking_generation(lib) -> None:
+def test_stacking_generation(lib: Library) -> None:
     """
     Balance on one node with one fixed demand and 2 generations with different costs, on 1 timestep.
     """
@@ -198,7 +199,7 @@ def test_stacking_generation(lib) -> None:
     assert problem.solver.Objective().Value() == 30 * 100 + 50 * 50
 
 
-def test_spillage(lib) -> None:
+def test_spillage(lib: Library) -> None:
     """
     Balance on one node with one fixed demand and 1 generation higher than demand and 1 timestep .
     """
@@ -238,7 +239,7 @@ def test_spillage(lib) -> None:
     assert problem.solver.Objective().Value() == 30 * 200 + 50 * 10
 
 
-def test_min_up_down_times(lib) -> None:
+def test_min_up_down_times(lib: Library) -> None:
     """
     Model on 3 time steps with one thermal generation and one demand on a single node.
         - Demand is the following time series : [500 MW, 0, 0]
@@ -332,7 +333,7 @@ def test_min_up_down_times(lib) -> None:
     assert problem.solver.Objective().Value() == pytest.approx(72000, abs=0.01)
 
 
-def test_changing_demand(lib) -> None:
+def test_changing_demand(lib: Library) -> None:
     """
     Model on 3 time steps simple production, demand
         - P_max = 500 MW
@@ -388,7 +389,7 @@ def test_changing_demand(lib) -> None:
     assert problem.solver.Objective().Value() == 40000
 
 
-def test_min_up_down_times_2(lib) -> None:
+def test_min_up_down_times_2(lib: Library) -> None:
     """
     Model on 3 time steps with one thermal generation and one demand on a single node.
         - Demand is the following time series : [500 MW, 0, 0]

tests/functional/test_performance.py

@@ -139,7 +139,7 @@ def test_large_sum_inside_model_with_sum_operator() -> None:
             float_variable(
                 "var",
                 lower_bound=literal(1),
-                upper_bound=literal(1),
+                upper_bound=literal(2),
                 structure=IndexingStructure(True, False),
             ),
         ],

tests/functional/test_xpansion.py

@@ -388,8 +388,14 @@ def test_generation_xpansion_two_time_steps_two_scenarios(
 
     output = OutputValues(problem)
     expected_output = OutputValues()
-    expected_output.component("G1").var("generation").value = [[0, 200], [0, 100]]
-    expected_output.component("CAND").var("generation").value = [[300, 300], [200, 300]]
+    expected_output.component("G1").var("generation").value = [
+        [0.0, 200.0],
+        [0.0, 100.0],
+    ]
+    expected_output.component("CAND").var("generation").value = [
+        [300.0, 300.0],
+        [200.0, 300.0],
+    ]
     expected_output.component("CAND").var("p_max").value = 300.0
 
     assert output == expected_output, f"Output differs from expected: {output}"