diff --git a/src/andromede/expression/context_adder.py b/src/andromede/expression/context_adder.py index 812e95f..07ee97f 100644 --- a/src/andromede/expression/context_adder.py +++ b/src/andromede/expression/context_adder.py @@ -18,6 +18,8 @@ ComponentVariableNode, ExpressionNode, ParameterNode, + ProblemParameterNode, + ProblemVariableNode, VariableNode, ) from .visitor import visit @@ -48,6 +50,16 @@ def comp_parameter(self, node: ComponentParameterNode) -> ExpressionNode: "This expression has already been associated to another component." ) + def pb_variable(self, node: ProblemVariableNode) -> ExpressionNode: + raise ValueError( + "This expression has already been associated to another component." + ) + + def pb_parameter(self, node: ProblemParameterNode) -> ExpressionNode: + raise ValueError( + "This expression has already been associated to another component." + ) + def add_component_context(id: str, expression: ExpressionNode) -> ExpressionNode: return visit(expression, ContextAdder(id)) diff --git a/src/andromede/expression/copy.py b/src/andromede/expression/copy.py index be61c1a..184a435 100644 --- a/src/andromede/expression/copy.py +++ b/src/andromede/expression/copy.py @@ -23,6 +23,8 @@ ParameterNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, ScenarioOperatorNode, TimeEvalNode, TimeShiftNode, @@ -58,6 +60,16 @@ def comp_variable(self, node: ComponentVariableNode) -> ExpressionNode: def comp_parameter(self, node: ComponentParameterNode) -> ExpressionNode: return ComponentParameterNode(node.component_id, node.name) + def pb_variable(self, node: ProblemVariableNode) -> ExpressionNode: + return ProblemVariableNode( + node.component_id, node.name, node.time_index, node.scenario_index + ) + + def pb_parameter(self, node: ProblemParameterNode) -> ExpressionNode: + return ProblemParameterNode( + node.component_id, node.name, node.time_index, node.scenario_index + ) + def time_shift(self, node: TimeShiftNode) -> ExpressionNode: return TimeShiftNode(visit(node.operand, self), visit(node.time_shift, self)) diff --git a/src/andromede/expression/degree.py b/src/andromede/expression/degree.py index 899042c..597d08e 100644 --- a/src/andromede/expression/degree.py +++ b/src/andromede/expression/degree.py @@ -17,6 +17,8 @@ ComponentVariableNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, TimeEvalNode, TimeShiftNode, TimeSumNode, @@ -77,6 +79,12 @@ def comp_variable(self, node: ComponentVariableNode) -> int: def comp_parameter(self, node: ComponentParameterNode) -> int: return 0 + def pb_variable(self, node: ProblemVariableNode) -> int: + return 1 + + def pb_parameter(self, node: ProblemParameterNode) -> int: + return 0 + def time_shift(self, node: TimeShiftNode) -> int: return visit(node.operand, self) diff --git a/src/andromede/expression/equality.py b/src/andromede/expression/equality.py index e26cf60..2998303 100644 --- a/src/andromede/expression/equality.py +++ b/src/andromede/expression/equality.py @@ -28,8 +28,12 @@ from andromede.expression.expression import ( AllTimeSumNode, BinaryOperatorNode, + ComponentParameterNode, + ComponentVariableNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, ScenarioOperatorNode, TimeEvalNode, TimeShiftNode, @@ -73,6 +77,22 @@ def visit(self, left: ExpressionNode, right: ExpressionNode) -> bool: return self.variable(left, right) if isinstance(left, ParameterNode) and isinstance(right, ParameterNode): return self.parameter(left, right) + if isinstance(left, ComponentVariableNode) and isinstance( + right, ComponentVariableNode + ): + return self.comp_variable(left, right) + if isinstance(left, ComponentParameterNode) and isinstance( + right, ComponentParameterNode + ): + return self.comp_parameter(left, right) + if isinstance(left, ProblemVariableNode) and isinstance( + right, ProblemVariableNode + ): + return self.problem_variable(left, right) + if isinstance(left, ProblemParameterNode) and isinstance( + right, ProblemParameterNode + ): + return self.problem_parameter(left, right) if isinstance(left, TimeShiftNode) and isinstance(right, TimeShiftNode): return self.time_shift(left, right) if isinstance(left, TimeEvalNode) and isinstance(right, TimeEvalNode): @@ -130,6 +150,36 @@ def variable(self, left: VariableNode, right: VariableNode) -> bool: def parameter(self, left: ParameterNode, right: ParameterNode) -> bool: return left.name == right.name + def comp_variable( + self, left: ComponentVariableNode, right: ComponentVariableNode + ) -> bool: + return left.name == right.name and left.component_id == right.component_id + + def comp_parameter( + self, left: ComponentParameterNode, right: ComponentParameterNode + ) -> bool: + return left.name == right.name and left.component_id == right.component_id + + def problem_variable( + self, left: ProblemVariableNode, right: ProblemVariableNode + ) -> bool: + return ( + left.name == right.name + and left.component_id == right.component_id + and left.time_index == right.time_index + and left.scenario_index == right.scenario_index + ) + + def problem_parameter( + self, left: ProblemParameterNode, right: ProblemParameterNode + ) -> bool: + return ( + left.name == right.name + and left.component_id == right.component_id + and left.time_index == right.time_index + and left.scenario_index == right.scenario_index + ) + def time_shift(self, left: TimeShiftNode, right: TimeShiftNode) -> bool: return self.visit(left.time_shift, right.time_shift) and self.visit( left.operand, right.operand diff --git a/src/andromede/expression/evaluate.py b/src/andromede/expression/evaluate.py index 642a539..a7f524e 100644 --- a/src/andromede/expression/evaluate.py +++ b/src/andromede/expression/evaluate.py @@ -20,6 +20,8 @@ ComponentVariableNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, TimeEvalNode, TimeShiftNode, TimeSumNode, @@ -33,6 +35,7 @@ ScenarioOperatorNode, VariableNode, ) +from .indexing import IndexingStructureProvider from .visitor import ExpressionVisitorOperations, visit @@ -58,11 +61,6 @@ def get_component_variable_value(self, component_id: str, name: str) -> float: def get_component_parameter_value(self, component_id: str, name: str) -> float: ... - # TODO: Should this really be an abstract method ? Or maybe, only the Provider in _make_value_provider should implement it. And the context attribute in the InstancesIndexVisitor is a ValueProvider that implements the parameter_is_constant_over_time method. Maybe create a child class of ValueProvider like TimeValueProvider ? - @abstractmethod - def parameter_is_constant_over_time(self, name: str) -> bool: - ... - @dataclass(frozen=True) class EvaluationContext(ValueProvider): @@ -86,9 +84,6 @@ def get_component_variable_value(self, component_id: str, name: str) -> float: def get_component_parameter_value(self, component_id: str, name: str) -> float: raise NotImplementedError() - def parameter_is_constant_over_time(self, name: str) -> bool: - raise NotImplementedError() - @dataclass(frozen=True) class EvaluationVisitor(ExpressionVisitorOperations[float]): @@ -117,6 +112,12 @@ def comp_parameter(self, node: ComponentParameterNode) -> float: def comp_variable(self, node: ComponentVariableNode) -> float: return self.context.get_component_variable_value(node.component_id, node.name) + def pb_parameter(self, node: ProblemParameterNode) -> float: + raise ValueError("Should not reach here.") + + def pb_variable(self, node: ProblemVariableNode) -> float: + raise ValueError("Should not reach here.") + def time_shift(self, node: TimeShiftNode) -> float: raise NotImplementedError() @@ -141,32 +142,3 @@ def port_field_aggregator(self, node: PortFieldAggregatorNode) -> float: def evaluate(expression: ExpressionNode, value_provider: ValueProvider) -> float: return visit(expression, EvaluationVisitor(value_provider)) - - -@dataclass(frozen=True) -class InstancesIndexVisitor(EvaluationVisitor): - """ - Evaluates an expression given as instances index which should have no variable and constant parameter values. - """ - - def variable(self, node: VariableNode) -> float: - raise ValueError("An instance index expression cannot contain variable") - - def parameter(self, node: ParameterNode) -> float: - if not self.context.parameter_is_constant_over_time(node.name): - raise ValueError( - "Parameter given in an instance index expression must be constant over time" - ) - return self.context.get_parameter_value(node.name) - - def time_shift(self, node: TimeShiftNode) -> float: - raise ValueError("An instance index expression cannot contain time shift") - - def time_eval(self, node: TimeEvalNode) -> float: - raise ValueError("An instance index expression cannot contain time eval") - - def time_sum(self, node: TimeSumNode) -> float: - raise ValueError("An instance index expression cannot contain time sum") - - def all_time_sum(self, node: AllTimeSumNode) -> float: - raise ValueError("An instance index expression cannot contain time sum") diff --git a/src/andromede/expression/evaluate_parameters.py b/src/andromede/expression/evaluate_parameters.py index f7c590d..e6e9cfe 100644 --- a/src/andromede/expression/evaluate_parameters.py +++ b/src/andromede/expression/evaluate_parameters.py @@ -12,9 +12,8 @@ from abc import ABC, abstractmethod from dataclasses import dataclass -from typing import List -from andromede.expression.evaluate import InstancesIndexVisitor, ValueProvider +from andromede.expression.evaluate import ValueProvider from .copy import CopyVisitor from .expression import ( @@ -59,39 +58,3 @@ def resolve_parameters( expression: ExpressionNode, parameter_provider: ParameterValueProvider ) -> ExpressionNode: return visit(expression, ParameterResolver(parameter_provider)) - - -def float_to_int(value: float) -> int: - if isinstance(value, int) or value.is_integer(): - return int(value) - else: - raise ValueError(f"{value} is not an integer.") - - -def evaluate_time_id(expr: ExpressionNode, value_provider: ValueProvider) -> int: - float_time_id = visit(expr, InstancesIndexVisitor(value_provider)) - try: - time_id = float_to_int(float_time_id) - except ValueError: - print(f"{expr} does not represent an integer time index.") - return time_id - - -# def get_time_ids_from_instances_index( -# instances_index: InstancesTimeIndex, value_provider: ValueProvider -# ) -> List[int]: -# time_ids = [] -# if isinstance(instances_index.expressions, list): # List[ExpressionNode] -# for expr in instances_index.expressions: -# time_ids.append(evaluate_time_id(expr, value_provider)) -# -# elif isinstance(instances_index.expressions, ExpressionRange): # ExpressionRange -# start_id = evaluate_time_id(instances_index.expressions.start, value_provider) -# stop_id = evaluate_time_id(instances_index.expressions.stop, value_provider) -# step_id = 1 -# if instances_index.expressions.step is not None: -# step_id = evaluate_time_id(instances_index.expressions.step, value_provider) -# # ExpressionRange includes stop_id whereas range excludes it -# time_ids = list(range(start_id, stop_id + 1, step_id)) -# -# return time_ids diff --git a/src/andromede/expression/expression.py b/src/andromede/expression/expression.py index 4ae42c1..84279ec 100644 --- a/src/andromede/expression/expression.py +++ b/src/andromede/expression/expression.py @@ -191,6 +191,102 @@ def comp_param(component_id: str, name: str) -> ComponentParameterNode: return ComponentParameterNode(component_id, name) +@dataclass(frozen=True) +class TimeIndex: + pass + + +@dataclass(frozen=True) +class NoTimeIndex(TimeIndex): + """ + Some values do not depend on the timestep, this index should be used for them + (think one variable for all timesteps). + """ + + pass + + +@dataclass(frozen=True) +class TimeShift(TimeIndex): + """ + Represents the current timestep + a shift. + + This should only be used for nodes that actually depend on the timestep, + never for time-independent nodes (constant parameters ...). + """ + + timeshift: int + + +@dataclass(frozen=True) +class TimeStep(TimeIndex): + """ + Represents a given timestep, independently of the current timestep. + + This should only be used for nodes that actually depend on the timestep, + never for time-independent nodes (constant parameters ...). + """ + + timestep: int + + +@dataclass(frozen=True) +class ScenarioIndex: + pass + + +@dataclass(frozen=True) +class NoScenarioIndex(ScenarioIndex): + """ + Some values do not depend on the scenario, this index should be used for them + (think one variable for all timesteps). + """ + + pass + + +@dataclass(frozen=True) +class CurrentScenarioIndex(ScenarioIndex): + """ + Represents the current scenario. + + This should only be used for nodes that actually depend on the scenario, + never for scenario-independent nodes (constant parameters ...). + """ + + pass + + +@dataclass(frozen=True) +class OneScenarioIndex(ScenarioIndex): + """ + Represents a given scenario out of all scenarios. + + This should only be used for nodes that actually depend on the scenario, + never for scenario-independent nodes (constant parameters ...). + """ + + scenario: int + + +@dataclass(frozen=True, eq=False) +class ProblemParameterNode(ExpressionNode): + """ + Represents one parameter of the optimization problem + """ + + component_id: str + name: str + time_index: TimeIndex + scenario_index: ScenarioIndex + + +def problem_param( + component_id: str, name: str, time_index: TimeIndex, scenario_index: ScenarioIndex +) -> ProblemParameterNode: + return ProblemParameterNode(component_id, name, time_index, scenario_index) + + @dataclass(frozen=True, eq=False) class ComponentVariableNode(ExpressionNode): """ @@ -209,6 +305,24 @@ def comp_var(component_id: str, name: str) -> ComponentVariableNode: return ComponentVariableNode(component_id, name) +@dataclass(frozen=True, eq=False) +class ProblemVariableNode(ExpressionNode): + """ + Represents one variable of the optimization problem + """ + + component_id: str + name: str + time_index: TimeIndex + scenario_index: ScenarioIndex + + +def problem_var( + component_id: str, name: str, time_index: TimeIndex, scenario_index: ScenarioIndex +) -> ProblemVariableNode: + return ProblemVariableNode(component_id, name, time_index, scenario_index) + + @dataclass(frozen=True, eq=False) class LiteralNode(ExpressionNode): value: float diff --git a/src/andromede/expression/indexing.py b/src/andromede/expression/indexing.py index d6635a3..4bc8471 100644 --- a/src/andromede/expression/indexing.py +++ b/src/andromede/expression/indexing.py @@ -12,6 +12,7 @@ from abc import ABC, abstractmethod from dataclasses import dataclass +from typing import List from andromede.expression.indexing_structure import IndexingStructure @@ -29,6 +30,8 @@ ParameterNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, ScenarioOperatorNode, TimeEvalNode, TimeShiftNode, @@ -74,21 +77,32 @@ def literal(self, node: LiteralNode) -> IndexingStructure: def negation(self, node: NegationNode) -> IndexingStructure: return visit(node.operand, self) - def addition(self, node: AdditionNode) -> IndexingStructure: - operands = [visit(o, self) for o in node.operands] - res = operands[0] - for o in node.operands[1:]: + def _combine(self, operands: List[ExpressionNode]) -> IndexingStructure: + if not operands: + return IndexingStructure(False, False) + res = visit(operands[0], self) + if res.is_time_scenario_varying(): + return res + for o in operands[1:]: res = res | visit(o, self) + if res.is_time_scenario_varying(): + return res return res + def addition(self, node: AdditionNode) -> IndexingStructure: + # performance note: + # here we don't need to visit all nodes, we can stop as soon as + # index is true/true + return self._combine(node.operands) + def multiplication(self, node: MultiplicationNode) -> IndexingStructure: - return visit(node.left, self) | visit(node.right, self) + return self._combine([node.left, node.right]) def division(self, node: DivisionNode) -> IndexingStructure: - return visit(node.left, self) | visit(node.right, self) + return self._combine([node.left, node.right]) def comparison(self, node: ComparisonNode) -> IndexingStructure: - return visit(node.left, self) | visit(node.right, self) + return self._combine([node.left, node.right]) def variable(self, node: VariableNode) -> IndexingStructure: time = self.context.get_variable_structure(node.name).time == True @@ -110,6 +124,16 @@ def comp_parameter(self, node: ComponentParameterNode) -> IndexingStructure: node.component_id, node.name ) + def pb_variable(self, node: ProblemVariableNode) -> IndexingStructure: + raise ValueError( + "Not relevant to compute indexation on already instantiated problem variables." + ) + + def pb_parameter(self, node: ProblemParameterNode) -> IndexingStructure: + raise ValueError( + "Not relevant to compute indexation on already instantiated problem parameters." + ) + def time_shift(self, node: TimeShiftNode) -> IndexingStructure: return visit(node.operand, self) diff --git a/src/andromede/expression/operators_expansion.py b/src/andromede/expression/operators_expansion.py new file mode 100644 index 0000000..13e2d40 --- /dev/null +++ b/src/andromede/expression/operators_expansion.py @@ -0,0 +1,240 @@ +import dataclasses +from dataclasses import dataclass +from typing import Callable, TypeVar, Union + +from andromede.expression import CopyVisitor, ExpressionNode, sum_expressions, visit +from andromede.expression.expression import ( + AllTimeSumNode, + ComponentParameterNode, + ComponentVariableNode, + CurrentScenarioIndex, + NoScenarioIndex, + NoTimeIndex, + OneScenarioIndex, + ProblemParameterNode, + ProblemVariableNode, + ScenarioOperatorNode, + TimeEvalNode, + TimeShift, + TimeShiftNode, + TimeStep, + TimeSumNode, + problem_param, + problem_var, +) +from andromede.expression.indexing import IndexingStructureProvider + +ExpressionEvaluator = Callable[[ExpressionNode], int] + + +@dataclass(frozen=True) +class ProblemDimensions: + """ + Dimensions for the simulation window + """ + + timesteps_count: int + scenarios_count: int + + +@dataclass(frozen=True) +class ProblemIndex: + """ + Index of an object in the simulation window. + """ + + timestep: int + scenario: int + + +@dataclass(frozen=True) +class OperatorsExpansion(CopyVisitor): + """ + Replaces aggregators (time sum, expectations ...) by their + arithmetic expansion. + + This will allow to easily translate it to a plain linear expression later on, + without complex handling of operators. + + The obtained expression only contains `ProblemVariableNode` for variables + and `ProblemParameterNode` parameters. + """ + + timesteps_count: int + scenarios_count: int + evaluator: ExpressionEvaluator + structure_provider: IndexingStructureProvider + + def comp_variable(self, node: ComponentVariableNode) -> ExpressionNode: + structure = self.structure_provider.get_component_variable_structure( + node.component_id, node.name + ) + time_index = TimeShift(0) if structure.time else NoTimeIndex() + scenario_index = ( + CurrentScenarioIndex() if structure.scenario else NoScenarioIndex() + ) + return problem_var(node.component_id, node.name, time_index, scenario_index) + + def comp_parameter(self, node: ComponentParameterNode) -> ExpressionNode: + structure = self.structure_provider.get_component_parameter_structure( + node.component_id, node.name + ) + time_index = TimeShift(0) if structure.time else NoTimeIndex() + scenario_index = ( + CurrentScenarioIndex() if structure.scenario else NoScenarioIndex() + ) + return problem_param(node.component_id, node.name, time_index, scenario_index) + + def time_shift(self, node: TimeShiftNode) -> ExpressionNode: + shift = self.evaluator(node.time_shift) + operand = visit(node.operand, self) + return apply_timeshift(operand, shift) + + def time_eval(self, node: TimeEvalNode) -> ExpressionNode: + timestep = self.evaluator(node.eval_time) + operand = visit(node.operand, self) + return apply_timestep(operand, timestep) + + def time_sum(self, node: TimeSumNode) -> ExpressionNode: + from_shift = self.evaluator(node.from_time) + to_shift = self.evaluator(node.to_time) + operand = visit(node.operand, self) + nodes = [] + for t in range(from_shift, to_shift + 1): + nodes.append(apply_timeshift(operand, t)) + return sum_expressions(nodes) + + def all_time_sum(self, node: AllTimeSumNode) -> ExpressionNode: + nodes = [] + operand = visit(node.operand, self) + for t in range(self.timesteps_count): + # if we sum previously "evaluated" variables for example x[0], it's ok + nodes.append(apply_timestep(operand, t, allow_existing=True)) + return sum_expressions(nodes) + + def scenario_operator(self, node: ScenarioOperatorNode) -> ExpressionNode: + if node.name != "Expectation": + raise ValueError(f"Scenario operator not supported: {node.name}") + nodes = [] + operand = visit(node.operand, self) + for t in range(self.scenarios_count): + nodes.append(apply_scenario(operand, t)) + return sum_expressions(nodes) / self.scenarios_count + + def pb_parameter(self, node: ProblemParameterNode) -> ExpressionNode: + raise ValueError("Should not reach") + + def pb_variable(self, node: ProblemVariableNode) -> ExpressionNode: + raise ValueError("Should not reach") + + +def expand_operators( + expression: ExpressionNode, + dimensions: ProblemDimensions, + evaluator: ExpressionEvaluator, + structure_provider: IndexingStructureProvider, +) -> ExpressionNode: + return visit( + expression, + OperatorsExpansion( + dimensions.timesteps_count, + dimensions.scenarios_count, + evaluator, + structure_provider, + ), + ) + + +TimeIndexedNode = TypeVar( + "TimeIndexedNode", bound=Union[ProblemParameterNode, ProblemVariableNode] +) + + +@dataclass(frozen=True) +class ApplyTimeShift(CopyVisitor): + """ + Shifts all underlying expressions. + """ + + timeshift: int + + def _apply_timeshift(self, node: TimeIndexedNode) -> TimeIndexedNode: + current_index = node.time_index + if isinstance(current_index, TimeShift): + return dataclasses.replace( + node, time_index=TimeShift(current_index.timeshift + self.timeshift) + ) + if isinstance(current_index, TimeStep): + return dataclasses.replace( + node, time_index=TimeStep(current_index.timestep + self.timeshift) + ) + if isinstance(current_index, NoTimeIndex): + return node + raise ValueError("Unknown time index type.") + + def pb_parameter(self, node: ProblemParameterNode) -> ProblemParameterNode: + return self._apply_timeshift(node) + + def pb_variable(self, node: ProblemVariableNode) -> ProblemVariableNode: + return self._apply_timeshift(node) + + +def apply_timeshift(expression: ExpressionNode, timeshift: int) -> ExpressionNode: + return visit(expression, ApplyTimeShift(timeshift)) + + +@dataclass(frozen=True) +class ApplyTimeStep(CopyVisitor): + """ + Applies timestep to all underlying expressions. + """ + + timestep: int + allow_existing: bool = False + + def _apply_timestep(self, node: TimeIndexedNode) -> TimeIndexedNode: + current_index = node.time_index + if isinstance(current_index, TimeShift): + return dataclasses.replace( + node, time_index=TimeStep(current_index.timeshift + self.timestep) + ) + if isinstance(current_index, TimeStep): + if not self.allow_existing: + raise ValueError( + "Cannot override a previously defined timestep (for example (x[0])[1])." + ) + return node + if isinstance(current_index, NoTimeIndex): + return node + raise ValueError("Unknown time index type.") + + def pb_parameter(self, node: ProblemParameterNode) -> ExpressionNode: + return self._apply_timestep(node) + + def pb_variable(self, node: ProblemVariableNode) -> ExpressionNode: + return self._apply_timestep(node) + + +def apply_timestep( + expression: ExpressionNode, timestep: int, allow_existing: bool = False +) -> ExpressionNode: + return visit(expression, ApplyTimeStep(timestep, allow_existing)) + + +@dataclass(frozen=True) +class ApplyScenario(CopyVisitor): + scenario: int + + def pb_parameter(self, node: ProblemParameterNode) -> ExpressionNode: + if isinstance(node.scenario_index, NoScenarioIndex): + return node + return dataclasses.replace(node, scenario_index=OneScenarioIndex(self.scenario)) + + def pb_variable(self, node: ProblemVariableNode) -> ExpressionNode: + if isinstance(node.scenario_index, NoScenarioIndex): + return node + return dataclasses.replace(node, scenario_index=OneScenarioIndex(self.scenario)) + + +def apply_scenario(expression: ExpressionNode, scenario: int) -> ExpressionNode: + return visit(expression, ApplyScenario(scenario)) diff --git a/src/andromede/expression/print.py b/src/andromede/expression/print.py index 807bb1e..469a24e 100644 --- a/src/andromede/expression/print.py +++ b/src/andromede/expression/print.py @@ -20,6 +20,8 @@ ExpressionNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, TimeEvalNode, TimeShiftNode, TimeSumNode, @@ -99,6 +101,14 @@ def comp_variable(self, node: ComponentVariableNode) -> str: def comp_parameter(self, node: ComponentParameterNode) -> str: return f"{node.component_id}.{node.name}" + def pb_variable(self, node: ProblemVariableNode) -> str: + # TODO + return f"{node.component_id}.{node.name}" + + def pb_parameter(self, node: ProblemParameterNode) -> str: + # TODO + return f"{node.component_id}.{node.name}" + def time_shift(self, node: TimeShiftNode) -> str: return f"({visit(node.operand, self)}.shift({visit(node.time_shift, self)}))" diff --git a/src/andromede/expression/visitor.py b/src/andromede/expression/visitor.py index fe7510d..351be9e 100644 --- a/src/andromede/expression/visitor.py +++ b/src/andromede/expression/visitor.py @@ -31,6 +31,8 @@ ParameterNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, ScenarioOperatorNode, TimeEvalNode, TimeShiftNode, @@ -90,6 +92,14 @@ def comp_parameter(self, node: ComponentParameterNode) -> T: def comp_variable(self, node: ComponentVariableNode) -> T: ... + @abstractmethod + def pb_parameter(self, node: ProblemParameterNode) -> T: + ... + + @abstractmethod + def pb_variable(self, node: ProblemVariableNode) -> T: + ... + @abstractmethod def time_shift(self, node: TimeShiftNode) -> T: ... @@ -135,6 +145,10 @@ def visit(root: ExpressionNode, visitor: ExpressionVisitor[T]) -> T: return visitor.comp_parameter(root) elif isinstance(root, ComponentVariableNode): return visitor.comp_variable(root) + elif isinstance(root, ProblemParameterNode): + return visitor.pb_parameter(root) + elif isinstance(root, ProblemVariableNode): + return visitor.pb_variable(root) elif isinstance(root, AdditionNode): return visitor.addition(root) elif isinstance(root, MultiplicationNode): diff --git a/src/andromede/model/model.py b/src/andromede/model/model.py index dc1ecc4..2ad59ce 100644 --- a/src/andromede/model/model.py +++ b/src/andromede/model/model.py @@ -39,6 +39,8 @@ ComponentVariableNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, ScenarioOperatorNode, TimeEvalNode, TimeShiftNode, @@ -268,6 +270,16 @@ def comp_variable(self, node: ComponentVariableNode) -> None: "Port definition must not contain a variable associated to a component." ) + def pb_parameter(self, node: ProblemParameterNode) -> None: + raise ValueError( + "Port definition must not contain a parameter associated to a component." + ) + + def pb_variable(self, node: ProblemVariableNode) -> None: + raise ValueError( + "Port definition must not contain a variable associated to a component." + ) + def time_shift(self, node: TimeShiftNode) -> None: visit(node.operand, self) diff --git a/src/andromede/simulation/linear_expression.py b/src/andromede/simulation/linear_expression.py index 134269d..4421761 100644 --- a/src/andromede/simulation/linear_expression.py +++ b/src/andromede/simulation/linear_expression.py @@ -14,11 +14,17 @@ Specific modelling for "instantiated" linear expressions, with only variables and literal coefficients. """ -from dataclasses import dataclass, field +import dataclasses +from dataclasses import dataclass from typing import Callable, Dict, List, Optional, TypeVar, Union -from andromede.expression.indexing_structure import IndexingStructure -from andromede.expression.scenario_operator import ScenarioOperator +from andromede.expression.expression import ( + OneScenarioIndex, + ScenarioIndex, + TimeIndex, + TimeShift, + TimeStep, +) T = TypeVar("T") @@ -119,8 +125,8 @@ class TermKey: component_id: str variable_name: str - time_expansion: TimeExpansion - scenario_operator: Optional[ScenarioOperator] + time_index: Optional[int] + scenario_index: Optional[int] def _str_for_coeff(coeff: float) -> str: @@ -132,6 +138,25 @@ def _str_for_coeff(coeff: float) -> str: return "{:+g}".format(coeff) +def _time_index_to_str(time_index: TimeIndex) -> str: + if isinstance(time_index, TimeShift): + if time_index.timeshift == 0: + return "t" + elif time_index.timeshift > 0: + return f"t + {time_index.timeshift}" + else: + return f"t - {-time_index.timeshift}" + if isinstance(time_index, TimeStep): + return f"{time_index.timestep}" + return "" + + +def _scenario_index_to_str(scenario_index: ScenarioIndex) -> str: + if isinstance(scenario_index, OneScenarioIndex): + return f"{scenario_index.scenario}" + return "" + + def _str_for_time_expansion(exp: TimeExpansion) -> str: if isinstance(exp, TimeShiftExpansion): return f".shift({exp.shift})" @@ -156,11 +181,8 @@ class Term: coefficient: float component_id: str variable_name: str - structure: IndexingStructure = field( - default=IndexingStructure(time=True, scenario=True) - ) - time_expansion: TimeExpansion = TimeExpansion() - scenario_operator: Optional[ScenarioOperator] = None + time_index: Optional[int] + scenario_index: Optional[int] # TODO: It may be useful to define __add__, __sub__, etc on terms, which should return a linear expression ? @@ -169,10 +191,14 @@ def is_zero(self) -> bool: def __str__(self) -> str: # Useful for debugging tests - result = _str_for_coeff(self.coefficient) + str(self.variable_name) - result += _str_for_time_expansion(self.time_expansion) - if self.scenario_operator is not None: - result += f".{str(self.scenario_operator)}" + return repr(self) + + def __repr__(self) -> str: + # Useful for debugging tests + result = ( + f"{_str_for_coeff(self.coefficient)}{self.component_id}.{self.variable_name}" + f"[{self.time_index},{self.scenario_index}]" + ) return result @@ -180,83 +206,44 @@ def generate_key(term: Term) -> TermKey: return TermKey( term.component_id, term.variable_name, - term.time_expansion, - term.scenario_operator, + term.time_index, + term.scenario_index, ) def _merge_dicts( lhs: Dict[TermKey, Term], rhs: Dict[TermKey, Term], - merge_func: Callable[[Term, Term], Term], - neutral: float, + merge_func: Callable[[Optional[Term], Optional[Term]], Term], ) -> Dict[TermKey, Term]: res = {} - for k, v in lhs.items(): - res[k] = merge_func( - v, - rhs.get( - k, - Term( - neutral, - v.component_id, - v.variable_name, - v.structure, - v.time_expansion, - v.scenario_operator, - ), - ), - ) - for k, v in rhs.items(): + for k, left in lhs.items(): + right = rhs.get(k, None) + res[k] = merge_func(left, right) + for k, right in rhs.items(): if k not in lhs: - res[k] = merge_func( - Term( - neutral, - v.component_id, - v.variable_name, - v.structure, - v.time_expansion, - v.scenario_operator, - ), - v, - ) + res[k] = merge_func(None, right) return res -def _merge_is_possible(lhs: Term, rhs: Term) -> None: - if lhs.component_id != rhs.component_id or lhs.variable_name != rhs.variable_name: - raise ValueError("Cannot merge terms for different variables") - if ( - lhs.time_expansion != rhs.time_expansion - or lhs.scenario_operator != rhs.scenario_operator - ): - raise ValueError("Cannot merge terms with different operators") - if lhs.structure != rhs.structure: - raise ValueError("Cannot merge terms with different structures") - - -def _add_terms(lhs: Term, rhs: Term) -> Term: - _merge_is_possible(lhs, rhs) - return Term( - lhs.coefficient + rhs.coefficient, - lhs.component_id, - lhs.variable_name, - lhs.structure, - lhs.time_expansion, - lhs.scenario_operator, - ) +def _add_terms(lhs: Optional[Term], rhs: Optional[Term]) -> Term: + if lhs is not None and rhs is not None: + return dataclasses.replace(rhs, coefficient=lhs.coefficient + rhs.coefficient) + elif lhs is not None and rhs is None: + return lhs + elif lhs is None and rhs is not None: + return rhs + raise ValueError("Cannot add 2 null terms.") -def _substract_terms(lhs: Term, rhs: Term) -> Term: - _merge_is_possible(lhs, rhs) - return Term( - lhs.coefficient - rhs.coefficient, - lhs.component_id, - lhs.variable_name, - lhs.structure, - lhs.time_expansion, - lhs.scenario_operator, - ) +def _substract_terms(lhs: Optional[Term], rhs: Optional[Term]) -> Term: + if lhs is not None and rhs is not None: + return dataclasses.replace(lhs, coefficient=lhs.coefficient - rhs.coefficient) + elif lhs is not None and rhs is None: + return lhs + elif lhs is None and rhs is not None: + return dataclasses.replace(rhs, coefficient=-rhs.coefficient) + raise ValueError("Cannot subtract 2 null terms.") class LinearExpression: @@ -316,14 +303,14 @@ def str_for_constant(self) -> str: else: return "{:+g}".format(self.constant) - def __str__(self) -> str: + def __repr__(self) -> str: # Useful for debugging tests result = "" if self.is_zero(): result += "0" else: for term in self.terms.values(): - result += str(term) + result += repr(term) result += self.str_for_constant() @@ -341,7 +328,7 @@ def __iadd__(self, rhs: "LinearExpression") -> "LinearExpression": if not isinstance(rhs, LinearExpression): return NotImplemented self.constant += rhs.constant - aggregated_terms = _merge_dicts(self.terms, rhs.terms, _add_terms, 0) + aggregated_terms = _merge_dicts(self.terms, rhs.terms, _add_terms) self.terms = aggregated_terms self.remove_zeros_from_terms() return self @@ -356,7 +343,7 @@ def __isub__(self, rhs: "LinearExpression") -> "LinearExpression": if not isinstance(rhs, LinearExpression): return NotImplemented self.constant -= rhs.constant - aggregated_terms = _merge_dicts(self.terms, rhs.terms, _substract_terms, 0) + aggregated_terms = _merge_dicts(self.terms, rhs.terms, _substract_terms) self.terms = aggregated_terms self.remove_zeros_from_terms() return self @@ -397,9 +384,8 @@ def __imul__(self, rhs: "LinearExpression") -> "LinearExpression": term.coefficient * const_expr.constant, term.component_id, term.variable_name, - term.structure, - term.time_expansion, - term.scenario_operator, + term.time_index, + term.scenario_index, ) _copy_expression(left_expr, self) return self @@ -428,9 +414,8 @@ def __itruediv__(self, rhs: "LinearExpression") -> "LinearExpression": term.coefficient / rhs.constant, term.component_id, term.variable_name, - term.structure, - term.time_expansion, - term.scenario_operator, + term.time_index, + term.scenario_index, ) return self diff --git a/src/andromede/simulation/linearize.py b/src/andromede/simulation/linearize.py index 8656590..045b8c5 100644 --- a/src/andromede/simulation/linearize.py +++ b/src/andromede/simulation/linearize.py @@ -10,138 +10,256 @@ # # This file is part of the Antares project. -import dataclasses +from abc import ABC, abstractmethod from dataclasses import dataclass -from typing import Optional +from typing import Any, Dict, List, Optional, Union -import andromede.expression.scenario_operator -from andromede.expression.evaluate import ValueProvider -from andromede.expression.evaluate_parameters import evaluate_time_id +from andromede.expression import ( + AdditionNode, + DivisionNode, + ExpressionVisitor, + MultiplicationNode, + NegationNode, +) from andromede.expression.expression import ( AllTimeSumNode, ComparisonNode, ComponentParameterNode, ComponentVariableNode, + CurrentScenarioIndex, ExpressionNode, LiteralNode, + NoScenarioIndex, + NoTimeIndex, + OneScenarioIndex, ParameterNode, PortFieldAggregatorNode, PortFieldNode, + ProblemParameterNode, + ProblemVariableNode, + ScenarioIndex, ScenarioOperatorNode, TimeEvalNode, + TimeIndex, + TimeShift, TimeShiftNode, + TimeStep, TimeSumNode, VariableNode, ) -from andromede.expression.indexing import IndexingStructureProvider -from andromede.expression.visitor import ExpressionVisitorOperations, T, visit -from andromede.simulation.linear_expression import ( - AllTimeExpansion, - LinearExpression, - Term, - TimeEvalExpansion, - TimeExpansion, - TimeShiftExpansion, - TimeSumExpansion, - generate_key, -) +from andromede.expression.visitor import visit +from andromede.simulation.linear_expression import LinearExpression, Term, TermKey -def _apply_time_expansion( - input: LinearExpression, time_expansion: TimeExpansion -) -> LinearExpression: - result_terms = {} - for term in input.terms.values(): - term_with_operator = dataclasses.replace( - term, time_expansion=term.time_expansion.apply(time_expansion) +class ParameterGetter(ABC): + @abstractmethod + def get_parameter_value( + self, + component_id: str, + parameter_name: str, + timestep: Optional[int], + scenario: Optional[int], + ) -> float: + pass + + +@dataclass +class MutableTerm: + coefficient: float + component_id: str + variable_name: str + time_index: Optional[int] + scenario_index: Optional[int] + + def to_key(self) -> TermKey: + return TermKey( + self.component_id, + self.variable_name, + self.time_index, + self.scenario_index, ) - result_terms[generate_key(term_with_operator)] = term_with_operator - # TODO: How can we apply a shift on a parameter ? It seems impossible for now as parameters must already be evaluated... - result_expr = LinearExpression(result_terms, input.constant) - return result_expr + def to_term(self) -> Term: + return Term( + self.coefficient, + self.component_id, + self.variable_name, + self.time_index, + self.scenario_index, + ) + + +@dataclass +class LinearExpressionData: + terms: List[MutableTerm] + constant: float + + def build(self) -> LinearExpression: + res_terms: Dict[TermKey, Any] = {} + for t in self.terms: + k = t.to_key() + if k in res_terms: + current_t = res_terms[k] + current_t.coefficient += t.coefficient + else: + res_terms[k] = t + for k, v in res_terms.items(): + res_terms[k] = v.to_term() + return LinearExpression(res_terms, self.constant) @dataclass(frozen=True) -class LinearExpressionBuilder(ExpressionVisitorOperations[LinearExpression]): +class LinearExpressionBuilder(ExpressionVisitor[LinearExpressionData]): """ Reduces a generic expression to a linear expression. - Parameters should have been evaluated first. + The input expression must respect the constraints of the output of + the operators expansion expression: + it must only contain `ProblemVariableNode` for variables + and `ProblemParameterNode` parameters. It cannot contain anymore + time aggregators or scenario aggregators, nor port-related nodes. """ - structure_provider: IndexingStructureProvider - value_provider: Optional[ValueProvider] = None + # TODO: linear expressions should be re-usable for different timesteps and scenarios + timestep: Optional[int] + scenario: Optional[int] + value_provider: Optional[ParameterGetter] = None + + def negation(self, node: NegationNode) -> LinearExpressionData: + operand = visit(node.operand, self) + operand.constant = -operand.constant + for t in operand.terms: + t.coefficient = -t.coefficient + return operand + + def addition(self, node: AdditionNode) -> LinearExpressionData: + operands = [visit(o, self) for o in node.operands] + terms = [] + constant: float = 0 + for o in operands: + constant += o.constant + terms.extend(o.terms) + return LinearExpressionData(terms=terms, constant=constant) + + def multiplication(self, node: MultiplicationNode) -> LinearExpressionData: + lhs = visit(node.left, self) + rhs = visit(node.right, self) + if not lhs.terms: + multiplier = lhs.constant + actual_expr = rhs + elif not rhs.terms: + multiplier = rhs.constant + actual_expr = lhs + else: + raise ValueError( + "At least one operand of a multiplication must be a constant expression." + ) + actual_expr.constant *= multiplier + for t in actual_expr.terms: + t.coefficient *= multiplier + return actual_expr + + def division(self, node: DivisionNode) -> LinearExpressionData: + lhs = visit(node.left, self) + rhs = visit(node.right, self) + if rhs.terms: + raise ValueError( + "The second operand of a division must be a constant expression." + ) + divider = rhs.constant + actual_expr = lhs + actual_expr.constant /= divider + for t in actual_expr.terms: + t.coefficient /= divider + return actual_expr + + def _get_timestep(self, time_index: TimeIndex) -> Optional[int]: + if isinstance(time_index, TimeShift): + if self.timestep is None: + raise ValueError("Cannot shift a time-independent expression.") + return self.timestep + time_index.timeshift + if isinstance(time_index, TimeStep): + return time_index.timestep + if isinstance(time_index, NoTimeIndex): + return None + else: + raise TypeError(f"Type {type(time_index)} is not a valid TimeIndex type.") + + def _get_scenario(self, scenario_index: ScenarioIndex) -> Optional[int]: + if isinstance(scenario_index, OneScenarioIndex): + return scenario_index.scenario + elif isinstance(scenario_index, CurrentScenarioIndex): + return self.scenario + elif isinstance(scenario_index, NoScenarioIndex): + return None + else: + raise TypeError( + f"Type {type(scenario_index)} is not a valid ScenarioIndex type." + ) - def literal(self, node: LiteralNode) -> LinearExpression: - return LinearExpression([], node.value) + def literal(self, node: LiteralNode) -> LinearExpressionData: + return LinearExpressionData([], node.value) - def comparison(self, node: ComparisonNode) -> LinearExpression: + def comparison(self, node: ComparisonNode) -> LinearExpressionData: raise ValueError("Linear expression cannot contain a comparison operator.") - def variable(self, node: VariableNode) -> LinearExpression: + def variable(self, node: VariableNode) -> LinearExpressionData: raise ValueError( "Variables need to be associated with their component ID before linearization." ) - def parameter(self, node: ParameterNode) -> LinearExpression: + def parameter(self, node: ParameterNode) -> LinearExpressionData: raise ValueError("Parameters must be evaluated before linearization.") - def comp_variable(self, node: ComponentVariableNode) -> LinearExpression: - return LinearExpression( + def comp_variable(self, node: ComponentVariableNode) -> LinearExpressionData: + raise ValueError( + "Variables need to be associated with their timestep/scenario before linearization." + ) + + def pb_variable(self, node: ProblemVariableNode) -> LinearExpressionData: + return LinearExpressionData( [ - Term( + MutableTerm( 1, node.component_id, node.name, - self.structure_provider.get_component_variable_structure( - node.component_id, node.name - ), + time_index=self._get_timestep(node.time_index), + scenario_index=self._get_scenario(node.scenario_index), ) ], 0, ) - def comp_parameter(self, node: ComponentParameterNode) -> LinearExpression: - raise ValueError("Parameters must be evaluated before linearization.") + def comp_parameter(self, node: ComponentParameterNode) -> LinearExpressionData: + raise ValueError( + "Parameters need to be associated with their timestep/scenario before linearization." + ) - def time_eval(self, node: TimeEvalNode) -> LinearExpression: - operand_expr = visit(node.operand, self) - eval_time = evaluate_time_id(node.eval_time, self._value_provider()) - time_expansion = TimeEvalExpansion(eval_time) - return _apply_time_expansion(operand_expr, time_expansion) - - def time_shift(self, node: TimeShiftNode) -> LinearExpression: - operand_expr = visit(node.operand, self) - time_shift = evaluate_time_id(node.time_shift, self._value_provider()) - time_expansion = TimeShiftExpansion(time_shift) - return _apply_time_expansion(operand_expr, time_expansion) - - def time_sum(self, node: TimeSumNode) -> LinearExpression: - operand_expr = visit(node.operand, self) - from_shift = evaluate_time_id(node.from_time, self._value_provider()) - to_shift = evaluate_time_id(node.to_time, self._value_provider()) - time_expansion = TimeSumExpansion(from_shift, to_shift) - if operand_expr.constant != 0: - # We could multiply by number of steps, but not very safe, it might depend on block bounds - # will be handled when refactoring for better parametrs handling - raise ValueError( - "Summing an expression containing a constant is not supported for now." - ) - return _apply_time_expansion(operand_expr, time_expansion) - - def all_time_sum(self, node: AllTimeSumNode) -> LinearExpression: - operand_expr = visit(node.operand, self) - time_expansion = AllTimeExpansion() - if operand_expr.constant != 0: - # We could multiply by number of steps if we had them, but we don't - # will be handled when refactoring for better parametrs handling - raise ValueError( - "Summing an expression containing a constant is not supported for now." - ) - return _apply_time_expansion(operand_expr, time_expansion) + def pb_parameter(self, node: ProblemParameterNode) -> LinearExpressionData: + # TODO SL: not the best place to do this. + # in the future, we should evaluate coefficients of variables as time vectors once for all timesteps + time_index = self._get_timestep(node.time_index) + scenario_index = self._get_scenario(node.scenario_index) + return LinearExpressionData( + [], + self._value_provider().get_parameter_value( + node.component_id, node.name, time_index, scenario_index + ), + ) - def _value_provider(self) -> ValueProvider: + def time_eval(self, node: TimeEvalNode) -> LinearExpressionData: + raise ValueError("Time operators need to be expanded before linearization.") + + def time_shift(self, node: TimeShiftNode) -> LinearExpressionData: + raise ValueError("Time operators need to be expanded before linearization.") + + def time_sum(self, node: TimeSumNode) -> LinearExpressionData: + raise ValueError("Time operators need to be expanded before linearization.") + + def all_time_sum(self, node: AllTimeSumNode) -> LinearExpressionData: + raise ValueError("Time operators need to be expanded before linearization.") + + def _value_provider(self) -> ParameterGetter: if self.value_provider is None: raise ValueError( "A value provider must be specified to linearize a time operator node." @@ -150,30 +268,15 @@ def _value_provider(self) -> ValueProvider: ) return self.value_provider - def scenario_operator(self, node: ScenarioOperatorNode) -> LinearExpression: - scenario_operator_cls = getattr( - andromede.expression.scenario_operator, node.name - ) - if scenario_operator_cls.degree() > 1: - raise ValueError( - f"Cannot linearize expression with a non-linear operator: {scenario_operator_cls.__name__}" - ) - - operand_expr = visit(node.operand, self) - result_terms = {} - for term in operand_expr.terms.values(): - term_with_operator = dataclasses.replace( - term, scenario_operator=scenario_operator_cls() - ) - result_terms[generate_key(term_with_operator)] = term_with_operator - - result_expr = LinearExpression(result_terms, operand_expr.constant) - return result_expr + def scenario_operator(self, node: ScenarioOperatorNode) -> LinearExpressionData: + raise ValueError("Scenario operators need to be expanded before linearization.") - def port_field(self, node: PortFieldNode) -> LinearExpression: + def port_field(self, node: PortFieldNode) -> LinearExpressionData: raise ValueError("Port fields must be replaced before linearization.") - def port_field_aggregator(self, node: PortFieldAggregatorNode) -> LinearExpression: + def port_field_aggregator( + self, node: PortFieldAggregatorNode + ) -> LinearExpressionData: raise ValueError( "Port fields aggregators must be replaced before linearization." ) @@ -181,9 +284,13 @@ def port_field_aggregator(self, node: PortFieldAggregatorNode) -> LinearExpressi def linearize_expression( expression: ExpressionNode, - structure_provider: IndexingStructureProvider, - value_provider: Optional[ValueProvider] = None, + timestep: Optional[int], + scenario: Optional[int], + value_provider: Optional[ParameterGetter] = None, ) -> LinearExpression: return visit( - expression, LinearExpressionBuilder(structure_provider, value_provider) - ) + expression, + LinearExpressionBuilder( + value_provider=value_provider, timestep=timestep, scenario=scenario + ), + ).build() diff --git a/src/andromede/simulation/optimization.py b/src/andromede/simulation/optimization.py index f72784b..9dce7c9 100644 --- a/src/andromede/simulation/optimization.py +++ b/src/andromede/simulation/optimization.py @@ -14,33 +14,25 @@ The optimization module contains the logic to translate the input model into a mathematical optimization problem. """ - +import itertools import math -from abc import ABC, abstractmethod from dataclasses import dataclass from enum import Enum from typing import Dict, Iterable, List, Optional import ortools.linear_solver.pywraplp as lp -from andromede.expression import ( - EvaluationVisitor, - ExpressionNode, - ParameterValueProvider, - ValueProvider, - resolve_parameters, - visit, -) +from andromede.expression import EvaluationVisitor, ExpressionNode, ValueProvider, visit from andromede.expression.context_adder import add_component_context from andromede.expression.indexing import IndexingStructureProvider, compute_indexation from andromede.expression.indexing_structure import IndexingStructure +from andromede.expression.operators_expansion import ProblemDimensions, expand_operators from andromede.expression.port_resolver import PortFieldKey, resolve_port -from andromede.expression.scenario_operator import Expectation from andromede.model.common import ValueType from andromede.model.constraint import Constraint from andromede.model.model import PortFieldId from andromede.simulation.linear_expression import LinearExpression, Term -from andromede.simulation.linearize import linearize_expression +from andromede.simulation.linearize import ParameterGetter, linearize_expression from andromede.simulation.strategy import MergedProblemStrategy, ModelSelectionStrategy from andromede.simulation.time_block import TimeBlock from andromede.study.data import DataBase @@ -62,8 +54,8 @@ class TimestepComponentVariableKey: def _get_parameter_value( context: "OptimizationContext", - block_timestep: int, - scenario: int, + block_timestep: Optional[int], + scenario: Optional[int], component_id: str, name: str, ) -> float: @@ -72,36 +64,10 @@ def _get_parameter_value( 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 - ) - - -class TimestepValueProvider(ABC): - """ - Interface which provides numerical values for individual timesteps. - """ - - @abstractmethod - def get_value(self, block_timestep: int, scenario: int) -> float: - raise NotImplementedError() - - def _make_value_provider( context: "OptimizationContext", - block_timestep: int, - scenario: int, - component: Component, + block_timestep: Optional[int], + scenario: Optional[int], ) -> ValueProvider: """ Create a value provider which takes its values from @@ -110,7 +76,7 @@ def _make_value_provider( Cannot evaluate expressions which contain variables. """ - class Provider(ValueProvider): + class Impl(ValueProvider): def get_component_variable_value(self, component_id: str, name: str) -> float: raise NotImplementedError( "Cannot provide variable value at problem build time." @@ -131,139 +97,18 @@ def get_parameter_value(self, name: str) -> float: "Parameter must be associated to its component before resolution." ) - def parameter_is_constant_over_time(self, name: str) -> bool: - return not component.model.parameters[name].structure.time - - return Provider() + return Impl() -@dataclass(frozen=True) -class ExpressionTimestepValueProvider(TimestepValueProvider): - context: "OptimizationContext" - component: Component - expression: ExpressionNode - - # OptimizationContext has knowledge of the block, so that get_value only needs block_timestep and scenario to get the correct data value - - def get_value(self, block_timestep: int, scenario: int) -> float: - param_value_provider = _make_value_provider( - self.context, block_timestep, scenario, self.component - ) - visitor = EvaluationVisitor(param_value_provider) - return visit(self.expression, visitor) - - -def _make_parameter_value_provider( +def _compute_expression_value( + expression: ExpressionNode, context: "OptimizationContext", - block_timestep: int, - scenario: int, -) -> ParameterValueProvider: - """ - A value provider which takes its values from - the parameter values as defined in the network data. - - Cannot evaluate expressions which contain variables. - """ - - class Provider(ParameterValueProvider): - def get_component_parameter_value(self, component_id: str, name: str) -> float: - return _get_parameter_value( - context, block_timestep, scenario, component_id, name - ) - - def get_parameter_value(self, name: str) -> float: - raise ValueError( - "Parameters should have been associated with their component before resolution." - ) - - return Provider() - - -def _make_data_structure_provider( - network: Network, component: Component -) -> IndexingStructureProvider: - """ - Retrieve information in data structure (parameter and variable) from the model - """ - - class Provider(IndexingStructureProvider): - def get_component_variable_structure( - self, component_id: str, name: str - ) -> IndexingStructure: - return network.get_component(component_id).model.variables[name].structure - - def get_component_parameter_structure( - self, component_id: str, name: str - ) -> IndexingStructure: - return network.get_component(component_id).model.parameters[name].structure - - def get_parameter_structure(self, name: str) -> IndexingStructure: - return component.model.parameters[name].structure - - def get_variable_structure(self, name: str) -> IndexingStructure: - return component.model.variables[name].structure - - return Provider() - - -@dataclass(frozen=True) -class ComponentContext: - """ - Helper class to fill the optimization problem with component-related equations and variables. - """ - - opt_context: "OptimizationContext" - component: Component - - def get_values(self, expression: ExpressionNode) -> TimestepValueProvider: - """ - The returned value provider will evaluate the provided expression. - """ - return ExpressionTimestepValueProvider( - self.opt_context, self.component, expression - ) - - def add_variable( - self, - block_timestep: int, - scenario: int, - model_var_name: str, - variable: lp.Variable, - ) -> None: - self.opt_context.register_component_variable( - block_timestep, scenario, self.component.id, model_var_name, variable - ) - - def get_variable( - self, block_timestep: int, scenario: int, variable_name: str - ) -> lp.Variable: - return self.opt_context.get_component_variable( - block_timestep, - scenario, - self.component.id, - variable_name, - self.component.model.variables[variable_name].structure, - ) - - def linearize_expression( - self, - block_timestep: int, - scenario: int, - expression: ExpressionNode, - ) -> LinearExpression: - parameters_valued_provider = _make_parameter_value_provider( - self.opt_context, block_timestep, scenario - ) - evaluated_expr = resolve_parameters(expression, parameters_valued_provider) - - value_provider = _make_value_provider( - self.opt_context, block_timestep, scenario, self.component - ) - structure_provider = _make_data_structure_provider( - self.opt_context.network, self.component - ) - - return linearize_expression(evaluated_expr, structure_provider, value_provider) + block_timestep: Optional[int], + scenario: Optional[int], +) -> float: + value_provider = _make_value_provider(context, block_timestep, scenario) + visitor = EvaluationVisitor(value_provider) + return visit(expression, visitor) class BlockBorderManagement(Enum): @@ -291,11 +136,21 @@ class SolverVariableInfo: is_in_objective: bool +def float_to_int(value: float) -> int: + if isinstance(value, int) or value.is_integer(): + return int(value) + else: + raise ValueError(f"{value} is not an integer.") + + class OptimizationContext: """ Helper class to build the optimization problem. - Maintains some mappings between model and solver objects. - Also provides navigation method in the model (components by node ...). + + - Maintains mappings between model and solver objects, + - Maintains mappings between port fields and expressions, + - Provides implementations of interfaces required by various visitors + used to transform expressions (values providers ...). """ def __init__( @@ -317,6 +172,10 @@ def __init__( PortFieldKey, List[ExpressionNode] ] = {} + self._constant_value_provider = self._make_constant_value_provider() + self._indexing_structure_provider = self._make_data_structure_provider() + self._parameter_getter = self._make_parameter_getter() + @property def network(self) -> Network: return self._network @@ -333,8 +192,21 @@ def connection_fields_expressions(self) -> Dict[PortFieldKey, List[ExpressionNod return self._connection_fields_expressions # TODO: Need to think about data processing when creating blocks with varying or inequal time steps length (aggregation, sum ?, mean of data ?) - def block_timestep_to_absolute_timestep(self, block_timestep: int) -> int: - return self._block.timesteps[block_timestep] + def block_timestep_to_absolute_timestep( + self, block_timestep: Optional[int] + ) -> Optional[int]: + """ + Timestep may be None for parameters or variables that don't depend on time. + """ + if block_timestep is None: + return None + return self._block.timesteps[self.get_actual_block_timestep(block_timestep)] + + def get_actual_block_timestep(self, block_timestep: int) -> int: + if self._border_management == BlockBorderManagement.CYCLE: + return block_timestep % self.block_length() + else: + raise NotImplementedError() @property def database(self) -> DataBase: @@ -352,23 +224,15 @@ def get_time_indices(self, index_structure: IndexingStructure) -> Iterable[int]: def get_scenario_indices(self, index_structure: IndexingStructure) -> Iterable[int]: return range(self.scenarios) if index_structure.scenario else range(1) - # TODO: API to improve, variable_structure guides which of the indices block_timestep and scenario should be used def get_component_variable( self, - block_timestep: int, - scenario: int, + block_timestep: Optional[int], + scenario: Optional[int], component_id: str, variable_name: str, - variable_structure: IndexingStructure, ) -> lp.Variable: - block_timestep = self._manage_border_timesteps(block_timestep) - - # TODO: Improve design, variable_structure defines indexing - if variable_structure.time == False: - block_timestep = 0 - if variable_structure.scenario == False: - scenario = 0 - + if block_timestep is not None: + block_timestep = self._manage_border_timesteps(block_timestep) return self._component_variables[ TimestepComponentVariableKey( component_id, variable_name, block_timestep, scenario @@ -382,8 +246,8 @@ def get_all_component_variables( def register_component_variable( self, - block_timestep: int, - scenario: int, + block_timestep: Optional[int], + scenario: Optional[int], component_id: str, variable_name: str, variable: lp.Variable, @@ -397,9 +261,6 @@ def register_component_variable( ) self._component_variables[key] = variable - def get_component_context(self, component: Component) -> ComponentContext: - return ComponentContext(self, component) - def register_connection_fields_expressions( self, component_id: str, @@ -412,25 +273,128 @@ def register_connection_fields_expressions( expression ) + def evaluate_time_bound(self, expression: ExpressionNode) -> int: + res = visit(expression, EvaluationVisitor(self._constant_value_provider)) + return float_to_int(res) -def _get_indexing( - constraint: Constraint, provider: IndexingStructureProvider -) -> IndexingStructure: - return ( - compute_indexation(constraint.expression, provider) - or compute_indexation(constraint.lower_bound, provider) - or compute_indexation(constraint.upper_bound, provider) - ) + def _make_data_structure_provider(self) -> IndexingStructureProvider: + """ + Retrieve information in data structure (parameter and variable) from the model + """ + network = self.network + + class Impl(IndexingStructureProvider): + def get_component_variable_structure( + self, component_id: str, name: str + ) -> IndexingStructure: + return ( + network.get_component(component_id).model.variables[name].structure + ) + + def get_component_parameter_structure( + self, component_id: str, name: str + ) -> IndexingStructure: + return ( + network.get_component(component_id).model.parameters[name].structure + ) + + def get_parameter_structure(self, name: str) -> IndexingStructure: + raise RuntimeError("Component context should have been initialized.") + + def get_variable_structure(self, name: str) -> IndexingStructure: + raise RuntimeError("Component context should have been initialized.") + + return Impl() + + def expand_operators(self, expression: ExpressionNode) -> ExpressionNode: + dimensions = ProblemDimensions(self.block_length(), self.scenarios) + time_bound_evaluator = self.evaluate_time_bound + return expand_operators( + expression, + dimensions, + time_bound_evaluator, + self._indexing_structure_provider, + ) + + def _make_parameter_getter(self) -> ParameterGetter: + ctxt = self + + class Impl(ParameterGetter): + def get_parameter_value( + self, + component_id: str, + parameter_name: str, + timestep: Optional[int], + scenario: Optional[int], + ) -> float: + return _get_parameter_value( + ctxt, + timestep, + scenario, + component_id, + parameter_name, + ) + + return Impl() + + def linearize_expression( + self, + expanded: ExpressionNode, + timestep: Optional[int] = None, + scenario: Optional[int] = None, + ) -> LinearExpression: + return linearize_expression( + expanded, timestep, scenario, self._parameter_getter + ) + + def compute_indexing(self, expression: ExpressionNode) -> IndexingStructure: + return compute_indexation(expression, self._indexing_structure_provider) + + def _make_constant_value_provider(self) -> ValueProvider: + """ + Value provider which only provides values for constant parameters + """ + context = self + network = self.network + + class Impl(ValueProvider): + def get_component_variable_value( + self, component_id: str, name: str + ) -> float: + raise NotImplementedError( + "Cannot provide variable value at problem build time." + ) + + def get_component_parameter_value( + self, component_id: str, name: str + ) -> float: + model = network.get_component(component_id).model + structure = model.parameters[name].structure + if structure.time or structure.scenario: + raise ValueError(f"Parameter {name} is not constant.") + return _get_parameter_value(context, None, None, component_id, name) + + def get_variable_value(self, name: str) -> float: + raise NotImplementedError( + "Cannot provide variable value at problem build time." + ) + + def get_parameter_value(self, name: str) -> float: + raise ValueError( + "Parameter must be associated to its component before resolution." + ) + + return Impl() -def _compute_indexing_structure( - context: ComponentContext, constraint: Constraint +def _compute_indexing( + context: OptimizationContext, constraint: Constraint ) -> IndexingStructure: - data_structure_provider = _make_data_structure_provider( - context.opt_context.network, context.component + return ( + context.compute_indexing(constraint.expression) + or context.compute_indexing(constraint.lower_bound) + or context.compute_indexing(constraint.upper_bound) ) - constraint_indexing = _get_indexing(constraint, data_structure_provider) - return constraint_indexing def _instantiate_model_expression( @@ -452,32 +416,41 @@ def _instantiate_model_expression( def _create_constraint( solver: lp.Solver, - context: ComponentContext, + context: OptimizationContext, constraint: Constraint, ) -> None: """ Adds a component-related constraint to the solver. """ - constraint_indexing = _compute_indexing_structure(context, constraint) - for block_timestep in context.opt_context.get_time_indices(constraint_indexing): - for scenario in context.opt_context.get_scenario_indices(constraint_indexing): + expanded = context.expand_operators(constraint.expression) + constraint_indexing = _compute_indexing(context, constraint) + + for block_timestep in context.get_time_indices(constraint_indexing): + for scenario in context.get_scenario_indices(constraint_indexing): linear_expr_at_t = context.linearize_expression( - block_timestep, scenario, constraint.expression + expanded, block_timestep, scenario ) + # What happens if there is some time_operator in the bounds ? constraint_data = ConstraintData( name=constraint.name, - lower_bound=context.get_values(constraint.lower_bound).get_value( - block_timestep, scenario + lower_bound=_compute_expression_value( + constraint.lower_bound, + context, + block_timestep, + scenario, ), - upper_bound=context.get_values(constraint.upper_bound).get_value( - block_timestep, scenario + upper_bound=_compute_expression_value( + constraint.upper_bound, + context, + block_timestep, + scenario, ), expression=linear_expr_at_t, ) make_constraint( solver, - context.opt_context, + context, block_timestep, scenario, constraint_data, @@ -488,39 +461,25 @@ def _create_objective( solver: lp.Solver, opt_context: OptimizationContext, component: Component, - component_context: ComponentContext, objective_contribution: ExpressionNode, ) -> None: instantiated_expr = _instantiate_model_expression( objective_contribution, component.id, opt_context ) - # We have already checked in the model creation that the objective contribution is neither indexed by time nor by scenario - linear_expr = component_context.linearize_expression(0, 0, instantiated_expr) + expanded = opt_context.expand_operators(instantiated_expr) + linear_expr = opt_context.linearize_expression(expanded) obj: lp.Objective = solver.Objective() for term in linear_expr.terms.values(): - # TODO : How to handle the scenario operator in a general manner ? - if isinstance(term.scenario_operator, Expectation): - weight = 1 / opt_context.scenarios - scenario_ids = range(opt_context.scenarios) - else: - weight = 1 - scenario_ids = range(1) - - for scenario in scenario_ids: - solver_vars = _get_solver_vars( - term, - opt_context, - 0, - scenario, - ) - - for solver_var in solver_vars: - opt_context._solver_variables[solver_var].is_in_objective = True - obj.SetCoefficient( - solver_var, - obj.GetCoefficient(solver_var) + weight * term.coefficient, - ) + solver_var = _get_solver_var( + term, + opt_context, + ) + opt_context._solver_variables[solver_var].is_in_objective = True + obj.SetCoefficient( + solver_var, + obj.GetCoefficient(solver_var) + term.coefficient, + ) # This should have no effect on the optimization obj.SetOffset(linear_expr.constant + obj.offset()) @@ -534,27 +493,16 @@ class ConstraintData: expression: LinearExpression -def _get_solver_vars( +def _get_solver_var( term: Term, context: OptimizationContext, - block_timestep: int, - scenario: int, -) -> List[lp.Variable]: - timesteps = term.time_expansion.get_timesteps( - block_timestep, context.block_length() +) -> lp.Variable: + return context.get_component_variable( + term.time_index, + term.scenario_index, + term.component_id, + term.variable_name, ) - solver_vars = [] - for t in timesteps: - solver_vars.append( - context.get_component_variable( - t, - scenario, - term.component_id, - term.variable_name, - term.structure, - ) - ) - return solver_vars def make_constraint( @@ -563,38 +511,27 @@ def make_constraint( block_timestep: int, scenario: int, data: ConstraintData, -) -> Dict[str, lp.Constraint]: +) -> None: """ Adds constraint to the solver. """ - solver_constraints = {} constraint_name = f"{data.name}_t{block_timestep}_s{scenario}" solver_constraint: lp.Constraint = solver.Constraint(constraint_name) constant: float = 0 for term in data.expression.terms.values(): - solver_vars = _get_solver_vars( + solver_var = _get_solver_var( term, context, - block_timestep, - scenario, ) - for solver_var in solver_vars: - coefficient = term.coefficient + solver_constraint.GetCoefficient( - solver_var - ) - solver_constraint.SetCoefficient(solver_var, coefficient) - # TODO: On pourrait aussi faire que l'objet Constraint n'ait pas de terme constant dans son expression et que les constantes soit déjà prises en compte dans les bornes, ça simplifierait le traitement ici - constant += data.expression.constant + coefficient = term.coefficient + solver_constraint.GetCoefficient(solver_var) + solver_constraint.SetCoefficient(solver_var, coefficient) + constant += data.expression.constant solver_constraint.SetBounds( data.lower_bound - constant, data.upper_bound - constant ) - # TODO: this dictionary does not make sense, we override the content when there are multiple instances - solver_constraints[constraint_name] = solver_constraint - return solver_constraints - class OptimizationProblem: name: str @@ -647,9 +584,20 @@ def _register_connection_fields_definitions(self) -> None: expression=instantiated_expression, ) + def _solver_variable_name( + self, component_id: str, var_name: str, t: Optional[int], s: Optional[int] + ) -> str: + scenario_suffix = f"_s{s}" if s is not None else "" + block_suffix = f"_t{t}" if t is not None 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 + return f"{component_id}_{var_name}{block_suffix}{scenario_suffix}" + def _create_variables(self) -> None: for component in self.context.network.all_components: - component_context = self.context.get_component_context(component) model = component.model for model_var in self.strategy.get_variables(model): @@ -666,72 +614,57 @@ def _create_variables(self) -> None: model_var.upper_bound, component.id, self.context ) - var_name: str = f"{model_var.name}" - component_prefix = f"{component.id}_" if component.id else "" + time_indices: Iterable[Optional[int]] = [None] + if var_indexing.time: + time_indices = self.context.get_time_indices(var_indexing) + scenario_indices: Iterable[Optional[int]] = [None] + if var_indexing.scenario: + scenario_indices = self.context.get_scenario_indices(var_indexing) + + for t, s in itertools.product(time_indices, scenario_indices): + lower_bound = -self.solver.infinity() + upper_bound = self.solver.infinity() + if instantiated_lb_expr: + lower_bound = _compute_expression_value( + instantiated_lb_expr, self.context, t, s + ) + if instantiated_ub_expr: + upper_bound = _compute_expression_value( + instantiated_ub_expr, self.context, t, s + ) - 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 "" + solver_var_name = self._solver_variable_name( + component.id, model_var.name, t, s ) - for scenario in self.context.get_scenario_indices(var_indexing): - lower_bound = -self.solver.infinity() - upper_bound = self.solver.infinity() - if instantiated_lb_expr: - lower_bound = component_context.get_values( - instantiated_lb_expr - ).get_value(block_timestep, scenario) - if instantiated_ub_expr: - upper_bound = component_context.get_values( - instantiated_ub_expr - ).get_value(block_timestep, scenario) - - scenario_suffix = ( - f"_s{scenario}" - if var_indexing.is_scenario_varying() - and (self.context.scenarios > 1) - else "" + 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}" ) - # 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 - 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, - 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 + 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, + solver_var_name, + ) + else: + solver_var = self.solver.NumVar( + lower_bound, + upper_bound, + solver_var_name, + ) + self.context.register_component_variable( + t, s, component.id, model_var.name, solver_var + ) def _create_constraints(self) -> None: for component in self.context.network.all_components: @@ -754,22 +687,19 @@ def _create_constraints(self) -> None: ) _create_constraint( self.solver, - self.context.get_component_context(component), + self.context, instantiated_constraint, ) def _create_objectives(self) -> None: for component in self.context.network.all_components: - component_context = self.context.get_component_context(component) model = component.model - for objective in self.strategy.get_objectives(model): if objective is not None: _create_objective( self.solver, self.context, component, - component_context, objective, ) diff --git a/src/andromede/simulation/output_values.py b/src/andromede/simulation/output_values.py index cf94e7c..c832148 100644 --- a/src/andromede/simulation/output_values.py +++ b/src/andromede/simulation/output_values.py @@ -121,7 +121,11 @@ def value(self, values: Union[float, List[float], List[List[float]]]) -> None: self._size = (size_s, size_t) - def _set(self, timestep: int, scenario: int, value: float) -> None: + def _set( + self, timestep: Optional[int], scenario: Optional[int], value: float + ) -> None: + timestep = 0 if timestep is None else timestep + scenario = 0 if scenario is None else scenario key = TimeScenarioIndex(timestep, scenario) if key not in self._value: size_s = max(self._size[0], scenario + 1) @@ -200,9 +204,6 @@ def _build_components(self) -> None: return for key, value in self.problem.context.get_all_component_variables().items(): - if (key.block_timestep is None) or (key.scenario is None): - continue - ( self.component(key.component_id) .var(str(key.variable_name)) diff --git a/src/andromede/study/data.py b/src/andromede/study/data.py index 9c011d7..e4b7547 100644 --- a/src/andromede/study/data.py +++ b/src/andromede/study/data.py @@ -38,8 +38,8 @@ class ScenarioIndex: @dataclass(frozen=True) class AbstractDataStructure(ABC): - def get_value(self, timestep: int, scenario: int) -> float: - return NotImplemented + def get_value(self, timestep: Optional[int], scenario: Optional[int]) -> float: + raise NotImplementedError() @abstractmethod def check_requirement(self, time: bool, scenario: bool) -> bool: @@ -54,7 +54,7 @@ def check_requirement(self, time: bool, scenario: bool) -> bool: class ConstantData(AbstractDataStructure): value: float - def get_value(self, timestep: int, scenario: int) -> float: + def get_value(self, timestep: Optional[int], scenario: Optional[int]) -> float: return self.value # ConstantData can be used for time varying or constant models @@ -74,7 +74,9 @@ class TimeSeriesData(AbstractDataStructure): time_series: Dict[TimeIndex, float] - def get_value(self, timestep: int, scenario: int) -> float: + def get_value(self, timestep: Optional[int], scenario: Optional[int]) -> float: + if timestep is None: + raise KeyError("Time series data requires a time index.") return self.time_series[TimeIndex(timestep)] def check_requirement(self, time: bool, scenario: bool) -> bool: @@ -94,7 +96,9 @@ class ScenarioSeriesData(AbstractDataStructure): scenario_series: Dict[ScenarioIndex, float] - def get_value(self, timestep: int, scenario: int) -> float: + def get_value(self, timestep: Optional[int], scenario: Optional[int]) -> float: + if scenario is None: + raise KeyError("Scenario series data requires a scenario index.") return self.scenario_series[ScenarioIndex(scenario)] def check_requirement(self, time: bool, scenario: bool) -> bool: @@ -143,7 +147,11 @@ class TimeScenarioSeriesData(AbstractDataStructure): time_scenario_series: pd.DataFrame scenarization: Optional[Scenarization] = None - def get_value(self, timestep: int, scenario: int) -> float: + def get_value(self, timestep: Optional[int], scenario: Optional[int]) -> float: + if timestep is None: + raise KeyError("Time scenario data requires a time index.") + if scenario is None: + raise KeyError("Time scenario data requires a scenario index.") if self.scenarization: scenario = self.scenarization.get_scenario_for_year(scenario) value = str(self.time_scenario_series.iloc[timestep, scenario]) diff --git a/tests/functional/test_andromede.py b/tests/functional/test_andromede.py index ec21cef..48f30c2 100644 --- a/tests/functional/test_andromede.py +++ b/tests/functional/test_andromede.py @@ -206,7 +206,7 @@ def test_variable_bound() -> None: 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", + match="Upper and lower bounds of variable G_generation_t0_s0 have the same value: 0", ): problem = build_problem(network, database, TimeBlock(1, [0]), 1) diff --git a/tests/functional/test_performance.py b/tests/functional/test_performance.py index 02c6ebb..53ffd85 100644 --- a/tests/functional/test_performance.py +++ b/tests/functional/test_performance.py @@ -9,7 +9,8 @@ # SPDX-License-Identifier: MPL-2.0 # # This file is part of the Antares project. - +import cProfile +from pstats import SortKey from typing import cast import pytest @@ -112,12 +113,13 @@ def test_large_sum_outside_model_with_loop() -> None: assert problem.solver.Objective().Value() == obj_coeff +# Takes 3 minutes with current implementation !! def test_large_sum_inside_model_with_sum_operator() -> None: """ Test performance when the problem involves an expression with a high number of terms. Here the objective function is the sum over nb_terms terms with the sum() operator inside the model """ - nb_terms = 10_000 + nb_terms = 3000 scenarios = 1 time_blocks = [TimeBlock(0, list(range(nb_terms)))] @@ -210,7 +212,7 @@ def test_basic_balance_on_whole_year() -> None: """ scenarios = 1 - horizon = 8760 + horizon = 10000 time_block = TimeBlock(1, list(range(horizon))) database = DataBase() @@ -233,7 +235,9 @@ def test_basic_balance_on_whole_year() -> None: network.connect(PortRef(demand, "balance_port"), PortRef(node, "balance_port")) network.connect(PortRef(gen, "balance_port"), PortRef(node, "balance_port")) - problem = build_problem(network, database, time_block, scenarios) + with cProfile.Profile() as pr: + problem = build_problem(network, database, time_block, scenarios) + pr.print_stats(sort=SortKey.CUMULATIVE) status = problem.solver.Solve() assert status == problem.solver.OPTIMAL diff --git a/tests/unittests/expressions/test_expressions.py b/tests/unittests/expressions/test_expressions.py index b852185..f8d644a 100644 --- a/tests/unittests/expressions/test_expressions.py +++ b/tests/unittests/expressions/test_expressions.py @@ -76,9 +76,6 @@ def get_component_variable_value(self, component_id: str, name: str) -> float: def get_component_parameter_value(self, component_id: str, name: str) -> float: return self.parameters[comp_key(component_id, name)] - def parameter_is_constant_over_time(self, name: str) -> bool: - raise NotImplementedError() - def test_comp_parameter() -> None: add_node = AdditionNode([LiteralNode(1), ComponentVariableNode("comp1", "x")]) diff --git a/tests/unittests/expressions/test_linear_expressions.py b/tests/unittests/expressions/test_linear_expressions.py index f948339..519f4ec 100644 --- a/tests/unittests/expressions/test_linear_expressions.py +++ b/tests/unittests/expressions/test_linear_expressions.py @@ -14,49 +14,15 @@ import pytest -from andromede.expression.expression import AllTimeSumNode -from andromede.expression.scenario_operator import Expectation -from andromede.simulation.linear_expression import ( - AllTimeExpansion, - LinearExpression, - Term, - TermKey, - TimeShiftExpansion, - TimeSumExpansion, -) +from andromede.simulation.linear_expression import LinearExpression, Term, TermKey @pytest.mark.parametrize( "term, expected", [ - (Term(1, "c", "x"), "+x"), - (Term(-1, "c", "x"), "-x"), - (Term(2.50, "c", "x"), "+2.5x"), - (Term(-3, "c", "x"), "-3x"), - ( - Term(-3, "c", "x", time_expansion=TimeShiftExpansion(-1)), - "-3x.shift(-1)", - ), - (Term(-3, "c", "x", time_expansion=AllTimeExpansion()), "-3x.sum()"), - ( - Term( - -3, - "c", - "x", - time_expansion=TimeSumExpansion(2, 3), - ), - "-3x.sum(2, 3)", - ), - (Term(-3, "c", "x", scenario_operator=Expectation()), "-3x.expec()"), ( - Term( - -3, - "c", - "x", - time_expansion=AllTimeExpansion(), - scenario_operator=Expectation(), - ), - "-3x.sum().expec()", + Term(-2.50, "c", "x", time_index=10, scenario_index=15), + "-2.5c.x[10,15]", ), ], ) @@ -68,16 +34,16 @@ def test_printing_term(term: Term, expected: str) -> None: "coeff, var_name, constant, expec_str", [ (0, "x", 0, "0"), - (1, "x", 0, "+x"), - (1, "x", 1, "+x+1"), - (3.7, "x", 1, "+3.7x+1"), + (1, "x", 0, "+c.x[0,0]"), + (1, "x", 1, "+c.x[0,0]+1"), + (3.7, "x", 1, "+3.7c.x[0,0]+1"), (0, "x", 1, "+1"), ], ) def test_affine_expression_printing_should_reflect_required_formatting( coeff: float, var_name: str, constant: float, expec_str: str ) -> None: - expr = LinearExpression([Term(coeff, "c", var_name)], constant) + expr = LinearExpression([Term(coeff, "c", var_name, 0, 0)], constant) assert str(expr) == expec_str @@ -97,8 +63,8 @@ def test_constant_expressions(lhs: LinearExpression, rhs: LinearExpression) -> N @pytest.mark.parametrize( "terms_dict, constant, exp_terms, exp_constant", [ - ({"x": Term(0, "c", "x")}, 1, {}, 1), - ({"x": Term(1, "c", "x")}, 1, {"x": Term(1, "c", "x")}, 1), + ({"x": Term(0, "c", "x", 0, 0)}, 1, {}, 1), + ({"x": Term(1, "c", "x", 0, 0)}, 1, {"x": Term(1, "c", "x", 0, 0)}, 1), ], ) def test_instantiate_linear_expression_from_dict( @@ -116,71 +82,29 @@ def test_instantiate_linear_expression_from_dict( "e1, e2, expected", [ ( - LinearExpression([Term(10, "c", "x")], 1), - LinearExpression([Term(5, "c", "x")], 2), - LinearExpression([Term(15, "c", "x")], 3), + LinearExpression([Term(10, "c", "x", 0, 0)], 1), + LinearExpression([Term(5, "c", "x", 0, 0)], 2), + LinearExpression([Term(15, "c", "x", 0, 0)], 3), ), ( - LinearExpression([Term(10, "c1", "x")], 1), - LinearExpression([Term(5, "c2", "x")], 2), - LinearExpression([Term(10, "c1", "x"), Term(5, "c2", "x")], 3), + LinearExpression([Term(10, "c1", "x", 0, 0)], 1), + LinearExpression([Term(5, "c2", "x", 0, 0)], 2), + LinearExpression([Term(10, "c1", "x", 0, 0), Term(5, "c2", "x", 0, 0)], 3), ), ( - LinearExpression([Term(10, "c", "x")], 0), - LinearExpression([Term(5, "c", "y")], 0), - LinearExpression([Term(10, "c", "x"), Term(5, "c", "y")], 0), + LinearExpression([Term(10, "c", "x", 0, 0)], 0), + LinearExpression([Term(5, "c", "y", 0, 0)], 0), + LinearExpression([Term(10, "c", "x", 0, 0), Term(5, "c", "y", 0, 0)], 0), ), ( - LinearExpression(), - LinearExpression( - [Term(10, "c", "x", time_expansion=TimeShiftExpansion(-1))] - ), - LinearExpression( - [Term(10, "c", "x", time_expansion=TimeShiftExpansion(-1))] - ), + LinearExpression([Term(10, "c", "x", 0, 0)], 0), + LinearExpression([Term(5, "c", "x", 1, 0)], 0), + LinearExpression([Term(10, "c", "x", 0, 0), Term(5, "c", "x", 1, 0)], 0), ), ( - LinearExpression(), - LinearExpression([Term(10, "c", "x", time_expansion=AllTimeExpansion())]), - LinearExpression([Term(10, "c", "x", time_expansion=AllTimeExpansion())]), - ), - ( - LinearExpression([Term(10, "c", "x")]), - LinearExpression( - [Term(10, "c", "x", time_expansion=TimeShiftExpansion(-1))] - ), - LinearExpression( - [ - Term(10, "c", "x"), - Term(10, "c", "x", time_expansion=TimeShiftExpansion(-1)), - ] - ), - ), - ( - LinearExpression([Term(10, "c", "x")]), - LinearExpression( - [ - Term( - 10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ) - ] - ), - LinearExpression( - [ - Term(10, "c", "x"), - Term( - 10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ), - ] - ), + LinearExpression([Term(10, "c", "x", 0, 0)], 0), + LinearExpression([Term(5, "c", "x", 0, 1)], 0), + LinearExpression([Term(10, "c", "x", 0, 0), Term(5, "c", "x", 0, 1)], 0), ), ], ) @@ -199,8 +123,8 @@ def test_addition_of_linear_expressions_with_different_number_of_instances_shoul def test_operation_that_leads_to_term_with_zero_coefficient_should_be_removed_from_terms() -> ( None ): - e1 = LinearExpression([Term(10, "c", "x")], 1) - e2 = LinearExpression([Term(10, "c", "x")], 2) + e1 = LinearExpression([Term(10, "c", "x", 0, 0)], 1) + e2 = LinearExpression([Term(10, "c", "x", 0, 0)], 2) e3 = e2 - e1 assert e3.terms == {} @@ -209,47 +133,20 @@ def test_operation_that_leads_to_term_with_zero_coefficient_should_be_removed_fr "e1, e2, expected", [ ( - LinearExpression([Term(10, "c", "x")], 3), + LinearExpression([Term(10, "c", "x", 1, 0)], 3), LinearExpression([], 2), - LinearExpression([Term(20, "c", "x")], 6), + LinearExpression([Term(20, "c", "x", 1, 0)], 6), ), ( - LinearExpression([Term(10, "c", "x")], 3), + LinearExpression([Term(10, "c", "x", 0, 1)], 3), LinearExpression([], 1), - LinearExpression([Term(10, "c", "x")], 3), + LinearExpression([Term(10, "c", "x", 0, 1)], 3), ), ( - LinearExpression([Term(10, "c", "x")], 3), + LinearExpression([Term(10, "c", "x", 0, 0)], 3), LinearExpression(), LinearExpression(), ), - ( - LinearExpression( - [ - Term( - 10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ) - ], - 3, - ), - LinearExpression([], 2), - LinearExpression( - [ - Term( - 20, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ) - ], - 6, - ), - ), ], ) def test_multiplication( @@ -260,8 +157,8 @@ def test_multiplication( def test_multiplication_of_two_non_constant_terms_should_raise_value_error() -> None: - e1 = LinearExpression([Term(10, "c", "x")], 0) - e2 = LinearExpression([Term(5, "c", "x")], 0) + e1 = LinearExpression([Term(10, "c", "x", 0, 0)], 0) + e2 = LinearExpression([Term(5, "c", "x", 0, 0)], 0) with pytest.raises(ValueError) as exc: _ = e1 * e2 assert str(exc.value) == "Cannot multiply two non constant expression" @@ -271,34 +168,8 @@ def test_multiplication_of_two_non_constant_terms_should_raise_value_error() -> "e1, expected", [ ( - LinearExpression([Term(10, "c", "x")], 5), - LinearExpression([Term(-10, "c", "x")], -5), - ), - ( - LinearExpression( - [ - Term( - 10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ) - ], - 5, - ), - LinearExpression( - [ - Term( - -10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ) - ], - -5, - ), + LinearExpression([Term(10, "c", "x", 1, 2)], 5), + LinearExpression([Term(-10, "c", "x", 1, 2)], -5), ), ], ) @@ -310,97 +181,31 @@ def test_negation(e1: LinearExpression, expected: LinearExpression) -> None: "e1, e2, expected", [ ( - LinearExpression([Term(10, "c", "x")], 1), - LinearExpression([Term(5, "c", "x")], 2), - LinearExpression([Term(5, "c", "x")], -1), + LinearExpression([Term(10, "c", "x", 0, 0)], 1), + LinearExpression([Term(5, "c", "x", 0, 0)], 2), + LinearExpression([Term(5, "c", "x", 0, 0)], -1), ), ( - LinearExpression([Term(10, "c1", "x")], 1), - LinearExpression([Term(5, "c2", "x")], 2), - LinearExpression([Term(10, "c1", "x"), Term(-5, "c2", "x")], -1), - ), - ( - LinearExpression([Term(10, "c", "x")], 0), - LinearExpression([Term(5, "c", "y")], 0), - LinearExpression([Term(10, "c", "x"), Term(-5, "c", "y")], 0), - ), - ( - LinearExpression(), + LinearExpression([Term(10, "c1", "x", 0, 0)], 1), + LinearExpression([Term(5, "c2", "x", 0, 0)], 2), LinearExpression( - [ - Term( - 10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - ) - ] - ), - LinearExpression( - [ - Term( - -10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - ) - ] + [Term(10, "c1", "x", 0, 0), Term(-5, "c2", "x", 0, 0)], -1 ), ), ( - LinearExpression(), - LinearExpression([Term(10, "c", "x", time_expansion=AllTimeExpansion())]), - LinearExpression([Term(-10, "c", "x", time_expansion=AllTimeExpansion())]), + LinearExpression([Term(10, "c", "x", 0, 0)], 0), + LinearExpression([Term(5, "c", "y", 0, 0)], 0), + LinearExpression([Term(10, "c", "x", 0, 0), Term(-5, "c", "y", 0, 0)], 0), ), ( - LinearExpression([Term(10, "c", "x")]), - LinearExpression( - [ - Term( - 10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - ) - ] - ), - LinearExpression( - [ - Term(10, "c", "x"), - Term( - -10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - ), - ] - ), + LinearExpression([Term(10, "c", "x", 0, 0)], 0), + LinearExpression([Term(5, "c", "x", 1, 0)], 0), + LinearExpression([Term(10, "c", "x", 0, 0), Term(-5, "c", "x", 1, 0)], 0), ), ( - LinearExpression([Term(10, "c", "x")]), - LinearExpression( - [ - Term( - 10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ) - ] - ), - LinearExpression( - [ - Term(10, "c", "x"), - Term( - -10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ), - ] - ), + LinearExpression([Term(10, "c", "x", 0, 0)], 0), + LinearExpression([Term(5, "c", "x", 0, 1)], 0), + LinearExpression([Term(10, "c", "x", 0, 0), Term(-5, "c", "x", 0, 1)], 0), ), ], ) @@ -414,41 +219,14 @@ def test_substraction( "e1, e2, expected", [ ( - LinearExpression([Term(10, "c", "x")], 15), + LinearExpression([Term(10, "c", "x", 0, 1)], 15), LinearExpression([], 5), - LinearExpression([Term(2, "c", "x")], 3), + LinearExpression([Term(2, "c", "x", 0, 1)], 3), ), ( - LinearExpression([Term(10, "c", "x")], 15), + LinearExpression([Term(10, "c", "x", 1, 0)], 15), LinearExpression([], 1), - LinearExpression([Term(10, "c", "x")], 15), - ), - ( - LinearExpression( - [ - Term( - 10, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ) - ], - 15, - ), - LinearExpression([], 5), - LinearExpression( - [ - Term( - 2, - "c", - "x", - time_expansion=TimeShiftExpansion(-1), - scenario_operator=Expectation(), - ) - ], - 3, - ), + LinearExpression([Term(10, "c", "x", 1, 0)], 15), ), ], ) @@ -459,7 +237,7 @@ def test_division( def test_division_by_zero_sould_raise_zero_division_error() -> None: - e1 = LinearExpression([Term(10, "c", "x")], 15) + e1 = LinearExpression([Term(10, "c", "x", 0, 0)], 15) e2 = LinearExpression() with pytest.raises(ZeroDivisionError) as exc: _ = e1 / e2 @@ -467,7 +245,7 @@ def test_division_by_zero_sould_raise_zero_division_error() -> None: def test_division_by_non_constant_expr_sould_raise_value_error() -> None: - e1 = LinearExpression([Term(10, "c", "x")], 15) + e1 = LinearExpression([Term(10, "c", "x", 0, 0)], 15) e2 = LinearExpression() with pytest.raises(ValueError) as exc: _ = e2 / e1 diff --git a/tests/unittests/expressions/test_linearization.py b/tests/unittests/expressions/test_linearization.py index 56b5aa9..e7626b6 100644 --- a/tests/unittests/expressions/test_linearization.py +++ b/tests/unittests/expressions/test_linearization.py @@ -1,105 +1,153 @@ +from unittest.mock import Mock + import pytest -from andromede.expression import ExpressionNode, param, var -from andromede.expression.expression import comp_var -from andromede.simulation.linear_expression import ( - AllTimeExpansion, - LinearExpression, - Term, +from andromede.expression import ExpressionNode, LiteralNode, literal, var +from andromede.expression.expression import ( + ComponentVariableNode, + CurrentScenarioIndex, + NoScenarioIndex, + TimeShift, + comp_param, + comp_var, + problem_var, +) +from andromede.expression.indexing import IndexingStructureProvider +from andromede.expression.indexing_structure import IndexingStructure +from andromede.expression.operators_expansion import ( + ProblemDimensions, + ProblemIndex, + expand_operators, ) -from andromede.simulation.linearize import linearize_expression +from andromede.simulation.linear_expression import LinearExpression, Term +from andromede.simulation.linearize import ParameterGetter, linearize_expression -from .test_expressions import ComponentEvaluationContext, StructureProvider +from .test_expressions import StructureProvider +P = comp_param("c", "p") +X = comp_var("c", "x") +Y = comp_var("c", "y") + + +def var_at(var: ComponentVariableNode, timestep, scenario) -> LinearExpression: + return LinearExpression( + terms=[ + Term( + 1, + var.component_id, + var.name, + time_index=timestep, + scenario_index=scenario, + ) + ], + constant=0, + ) -def test_linearization() -> None: - x = comp_var("c", "x") - expr = (5 * x + 3) / 2 - provider = StructureProvider() - assert linearize_expression(expr, provider) == LinearExpression( - [Term(2.5, "c", "x")], 1.5 - ) +def X_at(t: int = 0, s: int = 0) -> LinearExpression: + return var_at(X, timestep=t, scenario=s) - with pytest.raises(ValueError): - linearize_expression(param("p") * x, provider) +def Y_at(t: int = 0, s: int = 0) -> LinearExpression: + return var_at(Y, timestep=t, scenario=s) -def test_linearization_of_non_linear_expressions_should_raise_value_error() -> None: - x = var("x") - expr = x.variance() - provider = StructureProvider() - with pytest.raises(ValueError) as exc: - linearize_expression(expr, provider) - assert ( - str(exc.value) - == "Cannot linearize expression with a non-linear operator: Variance" - ) +def constant(c: float) -> LinearExpression: + return LinearExpression([], c) -def test_time_sum_is_distributed_on_expression() -> None: - x = comp_var("c", "x") - y = comp_var("c", "y") - expr = (x + y).time_sum() - provider = StructureProvider() +def evaluate_literal(node: ExpressionNode) -> int: + if isinstance(node, LiteralNode): + return int(node.value) + raise NotImplementedError("Can only evaluate literal nodes.") - assert linearize_expression(expr, provider) == LinearExpression( - [ - Term(1, "c", "x", time_expansion=AllTimeExpansion()), - Term(1, "c", "y", time_expansion=AllTimeExpansion()), - ], - 0, - ) +def test_linearization_before_operator_substitution_raises_an_error() -> None: + x = var("x") + expr = x.variance() -@pytest.mark.skip(reason="Not yet supported") -def test_time_sum_is_distributed_on_expression() -> None: - x = comp_var("c", "x") - y = comp_var("c", "y") - expr = (x + y).time_sum() provider = StructureProvider() + with pytest.raises( + ValueError, match="Scenario operators need to be expanded before linearization" + ): + linearize_expression(expr, timestep=0, scenario=0) - assert linearize_expression(expr, provider) == LinearExpression( - [ - Term(1, "c", "x", time_expansion=AllTimeExpansion()), - Term(1, "c", "y", time_expansion=AllTimeExpansion()), - ], - 0, - ) +class AllTimeScenarioDependent(IndexingStructureProvider): + def get_parameter_structure(self, name: str) -> IndexingStructure: + return IndexingStructure(True, True) -X = comp_var("c", "x") + def get_variable_structure(self, name: str) -> IndexingStructure: + return IndexingStructure(True, True) + def get_component_variable_structure( + self, component_id: str, name: str + ) -> IndexingStructure: + return IndexingStructure(True, True) -@pytest.mark.parametrize( - "expr", - [(X + 2).time_sum(), (X + 2).time_sum(-1, 2)], -) -def test_sum_of_constant_not_supported( + def get_component_parameter_structure( + self, component_id: str, name: str + ) -> IndexingStructure: + return IndexingStructure(True, True) + + +def _expand_and_linearize( expr: ExpressionNode, -) -> None: - structure_provider = StructureProvider() - value_provider = ComponentEvaluationContext() - with pytest.raises(ValueError): - linearize_expression(expr, structure_provider, value_provider) + dimensions: ProblemDimensions, + index: ProblemIndex, + parameter_value_provider: ParameterGetter, +) -> LinearExpression: + expanded = expand_operators( + expr, dimensions, evaluate_literal, AllTimeScenarioDependent() + ) + return linearize_expression( + expanded, index.timestep, index.scenario, parameter_value_provider + ) @pytest.mark.parametrize( - "expr", + "expr,expected", [ - X.shift(-1).shift(+1), - X.shift(-1).time_sum(), - X.shift(-1).time_sum(-2, +2), - X.time_sum().shift(-1), - X.time_sum(-2, +2).shift(-1), - X.eval(2).time_sum(), + ((5 * X + 3) / 2, constant(2.5) * X_at(t=0) + constant(1.5)), + ((X + Y).time_sum(), X_at(t=0) + Y_at(t=0) + X_at(t=1) + Y_at(t=1)), + (X.shift(-1).shift(+1), X_at(t=0)), + (X.shift(-1).time_sum(), X_at(t=-1) + X_at(t=0)), + (X.shift(-1).time_sum(-1, +1), X_at(t=-2) + X_at(t=-1) + X_at(t=0)), + (X.time_sum().shift(-1), X_at(t=-1) + X_at(t=0)), + (X.time_sum(-1, +1).shift(-1), X_at(t=-2) + X_at(t=-1) + X_at(t=0)), + (X.eval(2).time_sum(), X_at(t=2) + X_at(t=2)), + ((X + 2).time_sum(), X_at(t=0) + X_at(t=1) + constant(4)), + ((X + 2).time_sum(-1, 0), X_at(t=-1) + X_at(t=0) + constant(4)), + ((X + 2).time_sum(-1, 0), X_at(t=-1) + X_at(t=0) + constant(4)), ], ) -def test_linearization_of_nested_time_operations_should_raise_value_error( - expr: ExpressionNode, +def test_linearization_of_nested_time_operations( + expr: ExpressionNode, expected: LinearExpression ) -> None: - structure_provider = StructureProvider() - value_provider = ComponentEvaluationContext() - with pytest.raises(ValueError): - linearize_expression(expr, structure_provider, value_provider) + dimensions = ProblemDimensions(timesteps_count=2, scenarios_count=1) + index = ProblemIndex(timestep=0, scenario=0) + params = Mock(spec=ParameterGetter) + + assert _expand_and_linearize(expr, dimensions, index, params) == expected + + +def test_invalid_multiplication() -> None: + params = Mock(spec=ParameterGetter) + + x = problem_var( + "c", "x", time_index=TimeShift(0), scenario_index=CurrentScenarioIndex() + ) + expression = x * x + with pytest.raises(ValueError, match="constant"): + linearize_expression(expression, 0, 0, params) + + +def test_invalid_division() -> None: + params = Mock(spec=ParameterGetter) + + x = problem_var( + "c", "x", time_index=TimeShift(0), scenario_index=CurrentScenarioIndex() + ) + expression = literal(1) / x + with pytest.raises(ValueError, match="constant"): + linearize_expression(expression, 0, 0, params) diff --git a/tests/unittests/expressions/test_operators_expansion.py b/tests/unittests/expressions/test_operators_expansion.py new file mode 100644 index 0000000..05fff22 --- /dev/null +++ b/tests/unittests/expressions/test_operators_expansion.py @@ -0,0 +1,131 @@ +from dataclasses import dataclass +from typing import Dict + +import pytest + +from andromede.expression import ExpressionNode, LiteralNode +from andromede.expression.equality import expressions_equal +from andromede.expression.expression import ( + CurrentScenarioIndex, + NoScenarioIndex, + NoTimeIndex, + ProblemParameterNode, + ProblemVariableNode, + TimeShift, + TimeStep, + comp_param, + comp_var, + problem_param, + problem_var, +) +from andromede.expression.indexing import IndexingStructureProvider +from andromede.expression.indexing_structure import IndexingStructure +from andromede.expression.operators_expansion import ProblemDimensions, expand_operators + +P = comp_param("c", "p") +X = comp_var("c", "x") +CONST = comp_var("c", "const") + + +def shifted_P(t: int = 0) -> ProblemParameterNode: + return problem_param("c", "p", TimeShift(t), CurrentScenarioIndex()) + + +def P_at(t: int = 0) -> ProblemParameterNode: + return problem_param("c", "p", TimeStep(t), CurrentScenarioIndex()) + + +def X_at(t: int = 0) -> ProblemVariableNode: + return problem_var("c", "x", TimeStep(t), CurrentScenarioIndex()) + + +def shifted_X(t: int = 0) -> ProblemVariableNode: + return problem_var("c", "x", TimeShift(t), CurrentScenarioIndex()) + + +def const() -> ProblemVariableNode: + return problem_var("c", "x", NoTimeIndex(), NoScenarioIndex()) + + +def evaluate_literal(node: ExpressionNode) -> int: + if isinstance(node, LiteralNode): + return int(node.value) + raise NotImplementedError("Can only evaluate literal nodes.") + + +@dataclass(frozen=True) +class AllTimeScenarioDependent(IndexingStructureProvider): + def get_parameter_structure(self, name: str) -> IndexingStructure: + return IndexingStructure(True, True) + + def get_variable_structure(self, name: str) -> IndexingStructure: + return IndexingStructure(True, True) + + def get_component_variable_structure( + self, component_id: str, name: str + ) -> IndexingStructure: + return IndexingStructure(True, True) + + def get_component_parameter_structure( + self, component_id: str, name: str + ) -> IndexingStructure: + return IndexingStructure(True, True) + + +@dataclass(frozen=True) +class StructureProviderDict(IndexingStructureProvider): + """ + Defines indexing structure through dictionaries. Default is time-scenario dependent. + """ + + variables: Dict[str, IndexingStructure] + parameters: Dict[str, IndexingStructure] + + def get_parameter_structure(self, name: str) -> IndexingStructure: + return self.parameters.get(name, IndexingStructure(True, True)) + + def get_variable_structure(self, name: str) -> IndexingStructure: + return self.variables.get(name, IndexingStructure(True, True)) + + def get_component_variable_structure( + self, component_id: str, name: str + ) -> IndexingStructure: + return self.variables.get(name, IndexingStructure(True, True)) + + def get_component_parameter_structure( + self, component_id: str, name: str + ) -> IndexingStructure: + return self.parameters.get(name, IndexingStructure(True, True)) + + +@pytest.mark.parametrize( + "expr,expected", + [ + (X.time_sum(), X_at(0) + X_at(1)), + (X.shift(-1), shifted_X(-1)), + (X.time_sum(-2, 0), shifted_X(-2) + (shifted_X(-1) + shifted_X(0))), + ((P * X).shift(-1), shifted_P(-1) * shifted_X(-1)), + (X.shift(-1).shift(+1), shifted_X(0)), + ( + P * (P * X).time_sum(0, 1), + shifted_P(0) * (shifted_P(0) * shifted_X(0) + shifted_P(1) * shifted_X(1)), + ), + (X.eval(2).time_sum(), X_at(2) + X_at(2)), + ], +) +def test_operators_expansion(expr: ExpressionNode, expected: ExpressionNode) -> None: + expanded = expand_operators( + expr, ProblemDimensions(2, 1), evaluate_literal, AllTimeScenarioDependent() + ) + assert expressions_equal(expanded, expected) + + +def test_time_scenario_independent_var_has_no_time_or_scenario_index(): + structure_provider = StructureProviderDict( + parameters={}, variables={"const": IndexingStructure(False, False)} + ) + expr = (X + CONST).time_sum() + expanded = expand_operators( + expr, ProblemDimensions(2, 1), evaluate_literal, structure_provider + ) + assert expanded == X_at(0) + const() + X_at(1) + const()