From 4fc6b84fe3f17d6533765c5543192e773deb05a9 Mon Sep 17 00:00:00 2001 From: Bryn Pickering <17178478+brynpickering@users.noreply.github.com> Date: Fri, 19 Jul 2024 17:43:19 +0100 Subject: [PATCH] Feature: piecewise constraints (#569) * Only valid for two-dimensional SOS2 piecewise constraints. * `x` & `y` can be expressions (in pyomo) or pure decision variables (in pyomo & gurobi). --- CHANGELOG.md | 3 + docs/examples/piecewise_constraints.py | 271 ++++++++++++++++++ docs/migrating.md | 2 + docs/user_defined_math/components.md | 45 +++ .../examples/sos2_piecewise_linear_costs.yaml | 34 +++ mkdocs.yml | 1 + src/calliope/backend/backend_model.py | 229 ++++++++++++--- src/calliope/backend/expression_parser.py | 13 +- src/calliope/backend/gurobi_backend_model.py | 54 ++-- src/calliope/backend/latex_backend_model.py | 104 ++++++- src/calliope/backend/parsing.py | 74 +++-- src/calliope/backend/pyomo_backend_model.py | 133 ++++++--- src/calliope/backend/where_parser.py | 36 ++- src/calliope/config/math_schema.yaml | 24 ++ src/calliope/math/base.yaml | 4 +- .../sos2_piecewise_cost_investment.lp | 46 +++ tests/common/util.py | 3 +- tests/test_backend_general.py | 127 +++++++- tests/test_backend_gurobi.py | 90 ++++++ tests/test_backend_latex_backend.py | 84 +++++- tests/test_backend_pyomo.py | 116 ++++++++ tests/test_backend_where_parser.py | 5 +- tests/test_math.py | 35 +++ 23 files changed, 1346 insertions(+), 187 deletions(-) create mode 100644 docs/examples/piecewise_constraints.py create mode 100644 docs/user_defined_math/examples/sos2_piecewise_linear_costs.yaml create mode 100644 tests/common/lp_files/sos2_piecewise_cost_investment.lp diff --git a/CHANGELOG.md b/CHANGELOG.md index 4a1dc230..171b46c5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,9 @@ ### User-facing changes +|new| Piecewise constraints added to the YAML math with its own unique syntax (#107). +These constraints will be added to the optimisation problem using Special Ordered Sets of Type 2 (SOS2) variables. + |new| Direct interface to the Gurobi Python API using `!#yaml config.build.backend: gurobi` or `!#python model.build(backend="gurobi")`. Tests show that using the gurobi solver via the Python API reduces peak memory consumption and runtime by at least 30% for the combined model build and solve steps. This requires the `gurobipy` package which can be installed with `mamba`: `mamba install gurobi::gurobi`. diff --git a/docs/examples/piecewise_constraints.py b/docs/examples/piecewise_constraints.py new file mode 100644 index 00000000..ae6a44c5 --- /dev/null +++ b/docs/examples/piecewise_constraints.py @@ -0,0 +1,271 @@ +# --- +# jupyter: +# jupytext: +# custom_cell_magics: kql +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.16.1 +# kernelspec: +# display_name: calliope_docs_build [conda env:calliope-docs-new] +# language: python +# name: conda-env-calliope-docs-new-calliope_docs_build +# --- + +# %% [markdown] +# # Defining piecewise linear constraints +# +# In this tutorial, we use the national scale example model to implement a piecewise linear constraint. +# This constraint will represent a non-linear relationship between capacity and cost per unit capacity of Concentrating Solar Power (CSP). + +# %% + +import calliope +import numpy as np +import plotly.express as px + +calliope.set_log_verbosity("INFO", include_solver_output=False) + +# %% [markdown] +# # Model setup + +# %% [markdown] +# ## Defining our piecewise curve +# +# In the base national scale model, the CSP has a maximum rated capacity of 10,000 kW and a cost to invest in that capacity of 1000 USD / kW. +# +# In our updated model, the cost to invest in capacity will vary from 5000 USD / kW to 500 USD / kW as the CSP capacity increases: + +# %% +capacity_steps = [0, 2500, 5000, 7500, 10000] +cost_steps = [0, 3.75e6, 6e6, 7.5e6, 8e6] + +cost_per_cap = np.nan_to_num(np.divide(cost_steps, capacity_steps)).astype(int) + +fig = px.line( + x=capacity_steps, + y=cost_steps, + labels={"x": "Capacity (kW)", "y": "Investment cost (USD)"}, + markers="o", + range_y=[0, 10e6], + text=[f"{i} USD/kW" for i in cost_per_cap], +) +fig.update_traces(textposition="top center") +fig.show() + + +# %% [markdown] +# We can then provide this data when we load our model: +# +#
+#

Note

+#

+# We must index our piecewise data over "breakpoints". +#

+#
+# + +# %% +new_params = { + "parameters": { + "capacity_steps": { + "data": capacity_steps, + "index": [0, 1, 2, 3, 4], + "dims": "breakpoints", + }, + "cost_steps": { + "data": cost_steps, + "index": [0, 1, 2, 3, 4], + "dims": "breakpoints", + }, + } +} +print(new_params) +m = calliope.examples.national_scale(override_dict=new_params) + +# %% +m.inputs.capacity_steps + +# %% +m.inputs.cost_steps + +# %% [markdown] +# ## Creating our piecewise constraint +# +# We create the piecewise constraint by linking decision variables to the piecewise curve we have created. +# In this example, we require a new decision variable for investment costs that can take on the value defined by the curve at a given value of `flow_cap`. + +# %% +m.math["variables"]["piecewise_cost_investment"] = { + "description": "Investment cost that increases monotonically", + "foreach": ["nodes", "techs", "carriers", "costs"], + "where": "[csp] in techs", + "bounds": {"min": 0, "max": np.inf}, + "default": 0, +} + +# %% [markdown] +# We also need to link that decision variable to our total cost calculation. + +# %% +# Before +m.math["global_expressions"]["cost_investment_flow_cap"]["equations"] + +# %% +# Updated - we split the equation into two expressions. +m.math["global_expressions"]["cost_investment_flow_cap"]["equations"] = [ + {"expression": "$cost_sum * flow_cap", "where": "NOT [csp] in techs"}, + {"expression": "piecewise_cost_investment", "where": "[csp] in techs"}, +] + +# %% [markdown] +# We then need to define the piecewise constraint: + +# %% +m.math["piecewise_constraints"]["csp_piecewise_costs"] = { + "description": "Set investment costs values along a piecewise curve using special ordered sets of type 2 (SOS2).", + "foreach": ["nodes", "techs", "carriers", "costs"], + "where": "piecewise_cost_investment", + "x_expression": "flow_cap", + "x_values": "capacity_steps", + "y_expression": "piecewise_cost_investment", + "y_values": "cost_steps", +} + +# %% [markdown] +# Then we can build our optimisation problem: + +# %% [markdown] +# # Building and checking the optimisation problem +# +# With our piecewise constraint defined, we can build our optimisation problem + +# %% +m.build() + +# %% [markdown] +# And we can see that our piecewise constraint exists in the built optimisation problem "backend" + +# %% +m.backend.verbose_strings() +m.backend.get_piecewise_constraint("csp_piecewise_costs").to_series().dropna() + +# %% [markdown] +# ## Solve the optimisation problem +# +# Once we have all of our optimisation problem components set up as we desire, we can solve the problem. + +# %% +m.solve() + +# %% [markdown] +# The results are stored in `m._model_data` and can be accessed by the public property `m.results` + +# %% [markdown] +# ## Analysing the outputs + +# %% +# Absolute +csp_cost = m.results.cost_investment_flow_cap.sel(techs="csp") +csp_cost.to_series().dropna() + +# %% +# Relative to capacity +csp_cap = m.results.flow_cap.sel(techs="csp") +csp_cost_rel = csp_cost / csp_cap +csp_cost_rel.to_series().dropna() + +# %% +# Plotted on our piecewise curve +fig.add_scatter( + x=csp_cap.to_series().dropna().values, + y=csp_cost.to_series().dropna().values, + mode="markers", + marker_symbol="cross", + marker_size=10, + marker_color="cyan", + name="Installed capacity", +) +fig.show() + +# %% [markdown] +# ## YAML model definition +# We have updated the model parameters and math interactively in Python in this tutorial, the definition in YAML would look like: + +# %% [markdown] +# ### Math +# +# Saved as e.g., `csp_piecewise_math.yaml`. +# +# ```yaml +# variables: +# piecewise_cost_investment: +# description: Investment cost that increases monotonically +# foreach: [nodes, techs, carriers, costs] +# where: "[csp] in techs" +# bounds: +# min: 0 +# max: .inf +# default: 0 +# +# piecewise_constraints: +# csp_piecewise_costs: +# description: > +# Set investment costs values along a piecewise curve using special ordered sets of type 2 (SOS2). +# foreach: [nodes, techs, carriers, costs] +# where: "[csp] in techs" +# x_expression: flow_cap +# x_values: capacity_steps +# y_expression: piecewise_cost_investment +# y_values: cost_steps +# +# global_expressions: +# cost_investment_flow_cap.equations: +# - expression: "$cost_sum * flow_cap" +# where: "NOT [csp] in techs" +# - expression: "piecewise_cost_investment" +# where: "[csp] in techs" +# ``` + +# %% [markdown] +# ### Scenario definition +# +# Loaded into the national-scale example model with: `calliope.examples.national_scale(scenario="piecewise_csp_cost")` +# +# ```yaml +# overrides: +# piecewise_csp_cost: +# config.init.add_math: [csp_piecewise_math.yaml] +# parameters: +# capacity_steps: +# data: [0, 2500, 5000, 7500, 10000] +# index: [0, 1, 2, 3, 4] +# dims: "breakpoints" +# cost_steps: +# data: [0, 3.75e6, 6e6, 7.5e6, 8e6] +# index: [0, 1, 2, 3, 4] +# dims: "breakpoints" +# ``` + +# %% [markdown] +# ## Troubleshooting +# +# If you are failing to load a piecewise constraint or it isn't working as expected, here are some common things to note: +# +# 1. The extent of your `x_values` and `y_values` will dictate the maximum values of your piecewise decision variables. +# In this example, we define `capacity_steps` over the full capacity range that we allow our CSP to cover in the model. +# However, if we set `capacity_steps` to `[0, 2500, 5000, 7500, 9000]` then `flow_cap` would _never_ go above a value of 9000. +# +# 2. The `x_values` and `y_values` parameters must have the same number of breakpoints and be indexed over `breakpoints`. +# It is possible to extend these parameters to be indexed over other dimensions (e.g., different technologies with different piecewise curves) but it must _always_ include the `breakpoints` dimension. +# +# 3. `x_values` must increase monotonically. That is, `[0, 5000, 2500, 7500, 10000]` is not valid for `capacity_steps` in this example. +# `y_values`, on the other hand, _can_ vary any way you like; `[0, 6e6, 3.75e6, 8e6, 7.5e6]` is valid for `cost_steps`. +# +# 4. `x_expression` and `y_expression` _must_ include reference to at least one decision variable. +# It can be a math expression, not only a single decision variable. `flow_cap + storage_cap / 2` would be valid for `x_expression` in this example. +# +# 5. Piecewise constraints will make your problem more difficult to solve since each breakpoint adds a binary decision variable. +# Larger models with detailed piecewise constraints may not solve in a reasonable amount of time. +# diff --git a/docs/migrating.md b/docs/migrating.md index 6ce561c4..8dec46d2 100644 --- a/docs/migrating.md +++ b/docs/migrating.md @@ -977,5 +977,7 @@ Now, all components of our internal math are defined in a readable YAML syntax t You can add your own math to update the pre-defined math and to represent the physical system in ways we do not cover in our base math, or to apply new modelling methods and problem types (e.g., pathway or stochastic optimisation)! +When adding your own math, you can add [piecewise linear constraints](user_defined_math/components.md#piecewise-constraints), which is a new type of constraint compared to what could be defined in v0.6. + !!! info "See also" Our [pre-defined](pre_defined_math/index.md) and [user-defined](user_defined_math/index.md) math documentation. diff --git a/docs/user_defined_math/components.md b/docs/user_defined_math/components.md index f45e3e73..b180d347 100644 --- a/docs/user_defined_math/components.md +++ b/docs/user_defined_math/components.md @@ -87,6 +87,51 @@ Without a `where` string, all valid members (according to the `definition_matrix The equation expressions _must_ have comparison operators. 1. It can be deactivated so that it does not appear in the built optimisation problem by setting `active: false`. +## Piecewise constraints + +If you have non-linear relationships between two decision variables, you may want to represent them as a [piecewise linear function](https://en.wikipedia.org/wiki/Piecewise_linear_function). +The most common form of a piecewise function involves creating special ordered sets of type 2 (SOS2), set of binary variables that are linked together with specific constraints. + +!!! note + You can find a fully worked-out example in our [piecewise linear tutorial][defining-piecewise-linear-constraints]. + +Because the formulation of piecewise constraints is so specific, the math syntax differs from all other modelling components by having `x` and `y` attributes that need to be specified: + +```yaml +piecewise_constraints: + sos2_piecewise_flow_out: + description: Set outflow to follow capacity according to a piecewise curve. + foreach: [nodes, techs, carriers] + where: piecewise_x AND piecewise_y + x_expression: flow_cap + x_values: piecewise_x + y_expression: flow_out + y_values: piecewise_y + active: true +``` + +1. It needs a unique name (`sos2_piecewise_flow_out` in the above example). +1. Ideally, it has a long-form `description` added. +This is not required, but is useful metadata for later reference. +1. It can have a top-level `foreach` list and `where` string. +Without a `foreach`, it becomes an un-indexed constraint. +Without a `where` string, all valid members (according to the `definition_matrix`) based on `foreach` will be included in this constraint. +1. It has `x` and `y` [expression strings](syntax.md#expression-strings) (`x_expression`, `y_expression`). +1. It has `x` and `y` parameter references (`x_values`, `y_values`). +This should be a string name referencing an input parameter that contains the `breakpoints` dimension. +The values given by this parameter will be used to set the respective (`x` / `y`) expression at each breakpoint. +1. It can be deactivated so that it does not appear in the built optimisation problem by setting `active: false`. + +The component attributes combine to describe a piecewise curve that links the `x_expression` and `y_expression` according to their respective values in `x_values` and `y_values` at each breakpoint. + +!!! note + If the non-linear function you want to represent is convex, you may be able to avoid SOS2 variables, and instead represent it using [constraint components](#constraints). + You can find an example of this in our [piecewise linear costs custom math example][piecewise-linear-costs]. + +!!! warning + This approximation of a non-linear relationship may improve the representation of whatever real system you are modelling, but it will come at the cost of a more difficult model to solve. + Indeed, introducing piecewise constraints may mean your model can no longer reach a solution with the computational resources you have available. + ## Objectives With your constrained decision variables and a global expression that binds these variables to costs, you need an objective to minimise/maximise. The default, pre-defined objective is `min_cost_optimisation` and looks as follows: diff --git a/docs/user_defined_math/examples/sos2_piecewise_linear_costs.yaml b/docs/user_defined_math/examples/sos2_piecewise_linear_costs.yaml new file mode 100644 index 00000000..26e1edd0 --- /dev/null +++ b/docs/user_defined_math/examples/sos2_piecewise_linear_costs.yaml @@ -0,0 +1,34 @@ +# title: Piecewise linear costs - economies of scale +# +# description: | +# Add a piecewise cost function that reduces the incremental increase in investment costs with increasing technology rated capacity. +# This emulates "economies of scale", where the more of a technology there is deployed, the less expensive each additional investment in deployment. +# +# A more detailing example can be found in our [dedicated tutorial][defining-piecewise-linear-constraints]. +# +# New indexed parameters: +# +# - `piecewise_cost_investment_x` (defining the new set `breakpoints`) +# - `piecewise_cost_investment_y` (defining the new set `breakpoints`) +# +# --- + +variables: + piecewise_cost_investment: + description: Investment cost that increases monotonically + foreach: [nodes, techs, carriers, costs] + where: any(piecewise_cost_investment_x, over=breakpoints) AND any(piecewise_cost_investment_y, over=breakpoints) + bounds: + min: 0 + max: .inf + +piecewise_constraints: + sos2_piecewise_costs: + description: > + Set investment cost values along a piecewise curve using special ordered sets of type 2 (SOS2). + foreach: [nodes, techs, carriers, costs] + where: any(piecewise_cost_investment_x, over=breakpoints) AND any(piecewise_cost_investment_y, over=breakpoints) + x_expression: flow_cap + x_values: piecewise_cost_investment_x + y_expression: piecewise_cost_investment + y_values: piecewise_cost_investment_y diff --git a/mkdocs.yml b/mkdocs.yml index 5b64a556..b4e74cdd 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -128,6 +128,7 @@ nav: - examples/milp/index.md - examples/milp/notebook.py - examples/loading_tabular_data.py + - examples/piecewise_constraints.py - examples/calliope_model_object.py - examples/calliope_logging.py - Advanced features: diff --git a/src/calliope/backend/backend_model.py b/src/calliope/backend/backend_model.py index 40bff1e5..536938f8 100644 --- a/src/calliope/backend/backend_model.py +++ b/src/calliope/backend/backend_model.py @@ -25,7 +25,6 @@ ) import numpy as np -import pandas as pd import xarray as xr from calliope import exceptions @@ -48,7 +47,12 @@ T = TypeVar("T") _COMPONENTS_T = Literal[ - "variables", "constraints", "objectives", "parameters", "global_expressions" + "parameters", + "variables", + "global_expressions", + "constraints", + "piecewise_constraints", + "objectives", ] LOGGER = logging.getLogger(__name__) @@ -58,7 +62,14 @@ class BackendModelGenerator(ABC): """Helper class for backends.""" _VALID_COMPONENTS: tuple[_COMPONENTS_T, ...] = typing.get_args(_COMPONENTS_T) - _COMPONENT_ATTR_METADATA = ["description", "unit", "default", "title"] + _COMPONENT_ATTR_METADATA = [ + "description", + "unit", + "default", + "title", + "math_repr", + "original_dtype", + ] _PARAM_TITLES = extract_from_schema(MODEL_SCHEMA, "title") _PARAM_DESCRIPTIONS = extract_from_schema(MODEL_SCHEMA, "description") @@ -101,7 +112,7 @@ def add_parameter( @abstractmethod def add_constraint( - self, name: str, constraint_dict: parsing.UnparsedConstraintDict + self, name: str, constraint_dict: parsing.UnparsedConstraint ) -> None: """Add constraint equation to backend model in-place. @@ -110,13 +121,27 @@ def add_constraint( Args: name (str): Name of the constraint - constraint_dict (parsing.UnparsedConstraintDict): + constraint_dict (parsing.UnparsedConstraint): Constraint configuration dictionary, ready to be parsed and then evaluated. """ + @abstractmethod + def add_piecewise_constraint( + self, name: str, constraint_dict: parsing.UnparsedPiecewiseConstraint + ) -> None: + """Add piecewise constraint equation to backend model in-place. + + Resulting backend dataset entries will be piecewise constraint objects. + + Args: + name (str): + Name of the piecewise constraint + constraint_dict (parsing.UnparsedPiecewiseConstraint): Piecewise constraint configuration dictionary, ready to be parsed and then evaluated. + """ + @abstractmethod def add_global_expression( - self, name: str, expression_dict: parsing.UnparsedExpressionDict + self, name: str, expression_dict: parsing.UnparsedExpression ) -> None: """Add global expression (arithmetic combination of parameters and/or decision variables) to backend model in-place. @@ -124,25 +149,23 @@ def add_global_expression( Args: name (str): name of the global expression - expression_dict (parsing.UnparsedExpressionDict): global expression configuration dictionary, ready to be parsed and then evaluated. + expression_dict (parsing.UnparsedExpression): Global expression configuration dictionary, ready to be parsed and then evaluated. """ @abstractmethod - def add_variable( - self, name: str, variable_dict: parsing.UnparsedVariableDict - ) -> None: + def add_variable(self, name: str, variable_dict: parsing.UnparsedVariable) -> None: """Add decision variable to backend model in-place. Resulting backend dataset entries will be decision variable objects. Args: name (str): name of the variable. - variable_dict (parsing.UnparsedVariableDict): unparsed variable configuration dictionary. + variable_dict (parsing.UnparsedVariable): Unparsed variable configuration dictionary. """ @abstractmethod def add_objective( - self, name: str, objective_dict: parsing.UnparsedObjectiveDict + self, name: str, objective_dict: parsing.UnparsedObjective ) -> None: """Add objective arithmetic to backend model in-place. @@ -150,7 +173,7 @@ def add_objective( Args: name (str): name of the objective. - objective_dict (parsing.UnparsedObjectiveDict): unparsed objective configuration dictionary. + objective_dict (parsing.UnparsedObjective): Unparsed objective configuration dictionary. """ def log( @@ -207,12 +230,13 @@ def add_all_math(self): "variables", "global_expressions", "constraints", + "piecewise_constraints", "objectives", ]: component = components.removesuffix("s") - for name in self.inputs.math[components]: + for name, dict_ in self.inputs.math[components].items(): start = time.time() - getattr(self, f"add_{component}")(name) + getattr(self, f"add_{component}")(name, dict_) end = time.time() - start LOGGER.debug( f"Optimisation Model | {components}:{name} | Built in {end:.4f}s" @@ -242,10 +266,14 @@ def _add_run_mode_math(self) -> None: def _add_component( self, name: str, - component_dict: Tp | None, + component_dict: Tp, component_setter: Callable, component_type: Literal[ - "variables", "global_expressions", "constraints", "objectives" + "variables", + "global_expressions", + "constraints", + "piecewise_constraints", + "objectives", ], break_early: bool = True, ) -> parsing.ParsedBackendComponent | None: @@ -254,7 +282,7 @@ def _add_component( Args: name (str): name of the component. If not providing the `component_dict` directly, this name must be available in the input math provided on initialising the class. - component_dict (Tp | None): unparsed YAML dictionary configuration. + component_dict (Tp): unparsed YAML dictionary configuration. component_setter (Callable): function to combine evaluated xarray DataArrays into backend component objects. component_type (Literal["variables", "global_expressions", "constraints", "objectives"]): type of the added component. @@ -269,10 +297,8 @@ def _add_component( """ references: set[str] = set() - if component_dict is None: - component_dict = self.inputs.math[component_type][name] - if name not in self.inputs.math[component_type]: - self.inputs.math[component_type][name] = component_dict + if name not in self.inputs.math.get(component_type, {}): + self.inputs.math.set_key(f"{component_type}.name", component_dict) if break_early and not component_dict.get("active", True): self.log( @@ -424,30 +450,39 @@ def _add_to_dataset( All referenced objects will have their "references" attribute updated with this object's name. Defaults to None. """ - add_attrs = { - attr: unparsed_dict.pop(attr) - for attr in self._COMPONENT_ATTR_METADATA - if attr in unparsed_dict.keys() + yaml_snippet_attrs = {} + add_attrs = {} + for attr, val in unparsed_dict.items(): + if attr in self._COMPONENT_ATTR_METADATA: + add_attrs[attr] = val + else: + yaml_snippet_attrs[attr] = val + + if yaml_snippet_attrs: + add_attrs["yaml_snippet"] = AttrDict(yaml_snippet_attrs).to_yaml() + + da.attrs = { + "obj_type": obj_type, + "references": set(), + "coords_in_name": False, + **add_attrs, # type: ignore } - if unparsed_dict: - add_attrs["yaml_snippet"] = AttrDict(unparsed_dict).to_yaml() - - da.attrs.update( - { - "obj_type": obj_type, - "references": set(), - "coords_in_name": False, - **add_attrs, # type: ignore - } - ) self._dataset[name] = da - if references is not None: - for reference in references: - try: - self._dataset[reference].attrs["references"].add(name) - except KeyError: - continue + self._update_references(name, references) + + def _update_references(self, name: str, references: set): + """Update reference lists in dataset objects. + + Args: + name (str): Name to update in reference lists. + references (set): Names of dataset objects whose reference lists will be updated with `name`. + """ + for reference in references: + try: + self._dataset[reference].attrs["references"].add(name) + except KeyError: + continue @overload def _apply_func( # noqa: D102, override @@ -532,6 +567,11 @@ def constraints(self): """Slice of backend dataset to show only built constraints.""" return self._dataset.filter_by_attrs(obj_type="constraints") + @property + def piecewise_constraints(self): + """Slice of backend dataset to show only built piecewise constraints.""" + return self._dataset.filter_by_attrs(obj_type="piecewise_constraints") + @property def variables(self): """Slice of backend dataset to show only built variables.""" @@ -588,6 +628,79 @@ def __init__(self, inputs: xr.Dataset, instance: T, **kwargs) -> None: self.shadow_prices: ShadowPrices self._has_verbose_strings: bool = False + def add_piecewise_constraint( # noqa: D102, override + self, name: str, constraint_dict: parsing.UnparsedPiecewiseConstraint + ) -> None: + if "breakpoints" in constraint_dict.get("foreach", []): + raise BackendError( + f"(piecewise_constraints, {name}) | `breakpoints` dimension should not be in `foreach`. " + "Instead, index `x_values` and `y_values` parameters over `breakpoints`." + ) + + def _constraint_setter(where: xr.DataArray, references: set) -> xr.DataArray: + expressions = [] + vals = [] + for axis in ["x", "y"]: + expression_name = constraint_dict[f"{axis}_expression"] # type: ignore + parsed_component = parsing.ParsedBackendComponent( # type: ignore + "piecewise_constraints", + name, + {"equations": [{"expression": expression_name}], **constraint_dict}, + ) + eq = parsed_component.parse_equations(self.valid_component_names) + expression_da = eq[0].evaluate_expression( + self, where=where, references=references + ) + val_name = constraint_dict[f"{axis}_values"] # type: ignore + val_da = self.get_parameter(val_name) + if "breakpoints" not in val_da.dims: + raise BackendError( + f"(piecewise_constraints, {name}) | " + f"`{axis}_values` must be indexed over the `breakpoints` dimension." + ) + references.add(val_name) + expressions.append(expression_da) + vals.extend([*val_da.to_dataset("breakpoints").data_vars.values()]) + try: + return self._apply_func( + self._to_piecewise_constraint, + where, + 1, + *expressions, + *vals, + name=name, + n_breakpoints=len(self.inputs.breakpoints), + ) + except BackendError as err: + raise BackendError( + f"(piecewise_constraints, {name}) | Errors in generating piecewise constraint: {err}" + ) + + self._add_component( + name, constraint_dict, _constraint_setter, "piecewise_constraints" + ) + + @abstractmethod + def _to_piecewise_constraint( + self, x_var: Any, y_var: Any, *vals: float, name: str, n_breakpoints: int + ) -> Any: + """Utility function to generate a pyomo piecewise constraint for every element of an xarray DataArray. + + The x-axis decision variable need not be bounded. + This aligns piecewise constraint functionality with other possible backends (e.g., gurobipy). + + Args: + x_var (Any): The x-axis decision variable to constrain. + y_var (Any): The y-axis decision variable to constrain. + *vals (xr.DataArray): The x-axis and y-axis decision variable values at each piecewise constraint breakpoint. + name (str): The name of the piecewise constraint. + n_breakpoints (int): number of breakpoints + + Returns: + Any: + Return piecewise_constraint object. + """ + @abstractmethod def get_parameter(self, name: str, as_backend_objs: bool = True) -> xr.DataArray: """Extract parameter from backend dataset. @@ -642,6 +755,22 @@ def get_constraint( Otherwise, a xr.Dataset will be given, indexed over the same dimensions as the xr.DataArray, with variables for the constraint body, and upper (`ub`) and lower (`lb`) bounds. """ + def get_piecewise_constraint(self, name: str) -> xr.DataArray: + """Get piecewise constraint data as an array of backend interface objects. + + This method can be used to inspect and debug built piecewise constraints. + + Unlike other optimisation problem components, piecewise constraints can only be inspected as backend interface objects. + This is because each element is a collection of variables, parameters, constraints, and expressions. + + Args: + name (str): Name of piecewise constraint, as given in YAML piecewise constraint key. + + Returns: + xr.DataArray: Piecewise constraint array. + """ + return self._get_component(name, "piecewise_constraints") + @abstractmethod def get_variable(self, name: str, as_backend_objs: bool = True) -> xr.DataArray: """Extract decision variable array from backend dataset. @@ -909,7 +1038,15 @@ def _rebuild_references(self, references: set[str]) -> None: refs = [k for k in getattr(self, component).data_vars if k in references] for ref in refs: self.delete_component(ref, component) - getattr(self, "add_" + component.removesuffix("s"))(name=ref) + dict_ = self.inputs.attrs["math"][component][ref] + getattr(self, "add_" + component.removesuffix("s"))(ref, dict_) + + def _get_component(self, name: str, component_group: str) -> xr.DataArray: + component = getattr(self, component_group).get(name, None) + if component is None: + pretty_group_name = component_group.removesuffix("s").replace("_", " ") + raise KeyError(f"Unknown {pretty_group_name}: {name}") + return component def _get_variable_bound( self, bound: Any, name: str, references: set, fill_na: float | None = None @@ -958,9 +1095,7 @@ def _datetime_as_string(self, data: xr.DataArray | xr.Dataset) -> Iterator: yield finally: for name_ in datetime_coords: - data.coords[name_] = xr.apply_ufunc( - pd.to_datetime, data.coords[name_], keep_attrs=True - ) + data.coords[name_] = data.coords[name_].astype("datetime64[ns]") class ShadowPrices: diff --git a/src/calliope/backend/expression_parser.py b/src/calliope/backend/expression_parser.py index 6421e180..49139942 100644 --- a/src/calliope/backend/expression_parser.py +++ b/src/calliope/backend/expression_parser.py @@ -769,15 +769,12 @@ def as_math_string(self) -> str: # noqa: D102, override evaluated = self.as_array() self.eval_attrs["references"].add(self.name) - if evaluated.shape: - dims = rf"_\text{{{','.join(str(i).removesuffix('s') for i in evaluated.dims)}}}" + if evaluated.attrs["obj_type"] != "string": + data_var_string = evaluated.attrs["math_repr"] else: - dims = "" - if evaluated.attrs["obj_type"] in ["global_expressions", "variables"]: - formatted_name = rf"\textbf{{{self.name}}}" - elif evaluated.attrs["obj_type"] == "parameters": - formatted_name = rf"\textit{{{self.name}}}" - return formatted_name + dims + data_var_string = rf"\text{{{self.name}}}" + + return data_var_string def as_array(self) -> xr.DataArray: # noqa: D102, override backend_interface = self.eval_attrs["backend_interface"] diff --git a/src/calliope/backend/gurobi_backend_model.py b/src/calliope/backend/gurobi_backend_model.py index fab8ca27..e5f8096d 100644 --- a/src/calliope/backend/gurobi_backend_model.py +++ b/src/calliope/backend/gurobi_backend_model.py @@ -71,17 +71,17 @@ def add_parameter( # noqa: D102, override ) parameter_da = parameter_da.astype(float) - parameter_da.attrs["original_dtype"] = parameter_values.dtype attrs = { "title": self._PARAM_TITLES.get(parameter_name, None), "description": self._PARAM_DESCRIPTIONS.get(parameter_name, None), "unit": self._PARAM_UNITS.get(parameter_name, None), "default": default, + "original_dtype": parameter_values.dtype.name, } self._add_to_dataset(parameter_name, parameter_da, "parameters", attrs) def add_constraint( # noqa: D102, override - self, name: str, constraint_dict: parsing.UnparsedConstraintDict | None = None + self, name: str, constraint_dict: parsing.UnparsedConstraint ) -> None: def _constraint_setter( element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set @@ -94,7 +94,7 @@ def _constraint_setter( self._add_component(name, constraint_dict, _constraint_setter, "constraints") def add_global_expression( # noqa: D102, override - self, name: str, expression_dict: parsing.UnparsedExpressionDict | None = None + self, name: str, expression_dict: parsing.UnparsedExpression ) -> None: def _expression_setter( element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set @@ -110,11 +110,9 @@ def _expression_setter( ) def add_variable( # noqa: D102, override - self, name: str, variable_dict: parsing.UnparsedVariableDict | None = None + self, name: str, variable_dict: parsing.UnparsedVariable ) -> None: domain_dict = {"real": gurobipy.GRB.CONTINUOUS, "integer": gurobipy.GRB.INTEGER} - if variable_dict is None: - variable_dict = self.inputs.attrs["math"]["variables"][name] def _variable_setter(where: xr.DataArray, references: set): domain_type = domain_dict[variable_dict.get("domain", "real")] @@ -130,19 +128,15 @@ def _variable_setter(where: xr.DataArray, references: set): self._add_component(name, variable_dict, _variable_setter, "variables") def add_objective( # noqa: D102, override - self, name: str, objective_dict: parsing.UnparsedObjectiveDict | None = None + self, name: str, objective_dict: parsing.UnparsedObjective ) -> None: - min_ = gurobipy.GRB.MINIMIZE - max_ = gurobipy.GRB.MAXIMIZE sense_dict = { - "minimize": min_, - "minimise": min_, - "maximize": max_, - "maximise": max_, + "minimize": gurobipy.GRB.MINIMIZE, + "minimise": gurobipy.GRB.MINIMIZE, + "maximize": gurobipy.GRB.MAXIMIZE, + "maximise": gurobipy.GRB.MAXIMIZE, } - if objective_dict is None: - objective_dict = self.inputs.attrs["math"]["objectives"][name] sense = sense_dict[objective_dict["sense"]] def _objective_setter( @@ -285,7 +279,12 @@ def __renamer(val, *idx, name: str, attr: str): new_obj_name = f"{name}[{', '.join(idx)}]" setattr(val, attr, new_obj_name) - attribute_names = {"variables": "VarName", "constraints": "ConstrName"} + self._instance.update() + attribute_names = { + "variables": "VarName", + "constraints": "ConstrName", + "piecewise_constraints": "GenConstrName", + } with self._datetime_as_string(self._dataset): for da in self._dataset.filter_by_attrs(coords_in_name=False).values(): if da.attrs["obj_type"] not in attribute_names.keys(): @@ -300,7 +299,6 @@ def __renamer(val, *idx, name: str, attr: str): attr=attribute_names[da.attrs["obj_type"]], ) da.attrs["coords_in_name"] = True - self._instance.update() def to_lp(self, path: str | Path) -> None: # noqa: D102, override @@ -456,6 +454,28 @@ def has_integer_or_binary_variables(self) -> bool: # noqa: D102, override def _del_gurobi_obj(self, obj: Any) -> None: self._instance.remove(obj) + def _to_piecewise_constraint( # noqa: D102, override + self, + x_var: gurobipy.Var, + y_var: gurobipy.Var, + *vals: float, + name: str, + n_breakpoints: int, + ) -> gurobipy.GenConstr: + if not isinstance(x_var, gurobipy.Var) or not isinstance(y_var, gurobipy.Var): + raise BackendError( + "Gurobi backend can only build piecewise constraints using decision variables." + ) + y_vals = list(pd.Series(vals[n_breakpoints:]).dropna().values) + x_vals = list(pd.Series(vals[:n_breakpoints]).dropna().values) + try: + var = self._instance.addGenConstrPWL( + xpts=x_vals, ypts=y_vals, xvar=x_var, yvar=y_var, name=name + ) + except gurobipy.GurobiError as err: + raise BackendError(err) + return var + def _update_gurobi_variable( self, orig: gurobipy.Var, new: Any, *, bound: Literal["lb", "ub"] ) -> None: diff --git a/src/calliope/backend/latex_backend_model.py b/src/calliope/backend/latex_backend_model.py index 4e4d29d0..4dededdb 100644 --- a/src/calliope/backend/latex_backend_model.py +++ b/src/calliope/backend/latex_backend_model.py @@ -174,6 +174,10 @@ class LatexBackendModel(backend_model.BackendModelGenerator): Subject to ---------- + {% elif component_type == "piecewise_constraints" %} + + Subject to (piecewise) + ---------------------- {% elif component_type == "global_expressions" %} Where @@ -244,6 +248,8 @@ class LatexBackendModel(backend_model.BackendModelGenerator): \section{Objective} {% elif component_type == "constraints" %} \section{Subject to} + {% elif component_type == "piecewise_constraints" %} + \section{Subject to (piecewise)} {% elif component_type == "global_expressions" %} \section{Where} {% elif component_type == "variables" %} @@ -295,6 +301,9 @@ class LatexBackendModel(backend_model.BackendModelGenerator): {% elif component_type == "constraints" %} ## Subject to + {% elif component_type == "piecewise_constraints" %} + + ## Subject to (piecewise) {% elif component_type == "global_expressions" %} ## Where @@ -378,11 +387,14 @@ def add_parameter( # noqa: D102, override "title": self._PARAM_TITLES.get(parameter_name, None), "description": self._PARAM_DESCRIPTIONS.get(parameter_name, None), "unit": self._PARAM_UNITS.get(parameter_name, None), + "math_repr": rf"\textit{{{parameter_name}}}" + + self._dims_to_var_string(parameter_values), } + self._add_to_dataset(parameter_name, parameter_values, "parameters", attrs) def add_constraint( # noqa: D102, override - self, name: str, constraint_dict: parsing.UnparsedConstraintDict | None = None + self, name: str, constraint_dict: parsing.UnparsedConstraint | None = None ) -> None: equation_strings: list = [] @@ -400,8 +412,49 @@ def _constraint_setter( parsed_component, self.constraints[name], equations=equation_strings ) + def add_piecewise_constraint( # noqa: D102, override + self, name: str, constraint_dict: parsing.UnparsedPiecewiseConstraint + ) -> None: + non_where_refs: set = set() + + def _constraint_setter(where: xr.DataArray, references: set) -> xr.DataArray: + return where.where(where) + + math_parts = {} + for val in ["x_expression", "y_expression", "x_values", "y_values"]: + val_name = constraint_dict[val] + parsed_val = parsing.ParsedBackendComponent( + "piecewise_constraints", + name, + {"equations": [{"expression": val_name}]}, # type: ignore + ) + eq = parsed_val.parse_equations(self.valid_component_names) + math_parts[val] = eq[0].evaluate_expression( + self, return_type="math_string", references=non_where_refs + ) + + equation = { + "expression": rf"{ math_parts['y_expression']}\mathord{{=}}{ math_parts['y_values'] }", + "where": rf"{ math_parts['x_expression'] }\mathord{{=}}{math_parts['x_values']}", + } + if "foreach" in constraint_dict: + constraint_dict["foreach"].append("breakpoints") + else: + constraint_dict["foreach"] = ["breakpoints"] + parsed_component = self._add_component( + name, + constraint_dict, + _constraint_setter, + "piecewise_constraints", + break_early=False, + ) + self._update_references(name, non_where_refs) + self._generate_math_string( + parsed_component, self.piecewise_constraints[name], equations=[equation] + ) + def add_global_expression( # noqa: D102, override - self, name: str, expression_dict: parsing.UnparsedExpressionDict | None = None + self, name: str, expression_dict: parsing.UnparsedExpression | None = None ) -> None: equation_strings: list = [] @@ -418,36 +471,46 @@ def _expression_setter( "global_expressions", break_early=False, ) + expr_da = self.global_expressions[name] + expr_da.attrs["math_repr"] = rf"\textbf{{{name}}}" + self._dims_to_var_string( + expr_da + ) self._generate_math_string( - parsed_component, self.global_expressions[name], equations=equation_strings + parsed_component, expr_da, equations=equation_strings ) def add_variable( # noqa: D102, override - self, name: str, variable_dict: parsing.UnparsedVariableDict | None = None + self, name: str, variable_dict: parsing.UnparsedVariable | None = None ) -> None: domain_dict = {"real": r"\mathbb{R}\;", "integer": r"\mathbb{Z}\;"} + bound_refs: set = set() def _variable_setter(where: xr.DataArray, references: set) -> xr.DataArray: return where.where(where) - if variable_dict is None: - variable_dict = self.inputs.attrs["math"]["variables"][name] + domain = domain_dict[variable_dict.get("domain", "real")] parsed_component = self._add_component( name, variable_dict, _variable_setter, "variables", break_early=False ) - where_array = self.variables[name] + var_da = self.variables[name] + var_da.attrs["math_repr"] = rf"\textbf{{{name}}}" + self._dims_to_var_string( + var_da + ) domain = domain_dict[variable_dict.get("domain", "real")] - lb, ub = self._get_variable_bounds_string(name, variable_dict["bounds"]) + lb, ub = self._get_variable_bounds_string( + name, variable_dict["bounds"], bound_refs + ) + self._update_references(name, bound_refs.difference(name)) self._generate_math_string( - parsed_component, where_array, equations=[lb, ub], sense=r"\forall" + domain + parsed_component, var_da, equations=[lb, ub], sense=r"\forall" + domain ) def add_objective( # noqa: D102, override - self, name: str, objective_dict: parsing.UnparsedObjectiveDict | None = None + self, name: str, objective_dict: parsing.UnparsedObjective | None = None ) -> None: sense_dict = { "minimize": r"\min{}", @@ -455,8 +518,7 @@ def add_objective( # noqa: D102, override "minimise": r"\min{}", "maximise": r"\max{}", } - if objective_dict is None: - objective_dict = self.inputs.attrs["math"]["objectives"][name] + equation_strings: list = [] def _objective_setter( @@ -523,6 +585,7 @@ def generate_math_doc( for objtype in [ "objectives", "constraints", + "piecewise_constraints", "global_expressions", "variables", "parameters", @@ -607,10 +670,10 @@ def __mathify_text_in_text(instring): return jinja_env.from_string(template).render(**kwargs) def _get_variable_bounds_string( - self, name: str, bounds: parsing.UnparsedVariableBoundDict + self, name: str, bounds: parsing.UnparsedVariableBound, references: set ) -> tuple[dict[str, str], ...]: """Convert variable upper and lower bounds into math string expressions.""" - bound_dict: parsing.UnparsedConstraintDict = { + bound_dict: parsing.UnparsedConstraint = { "equations": [ {"expression": f"{bounds['min']} <= {name}"}, {"expression": f"{name} <= {bounds['max']}"}, @@ -619,6 +682,17 @@ def _get_variable_bounds_string( parsed_bounds = parsing.ParsedBackendComponent("constraints", name, bound_dict) equations = parsed_bounds.parse_equations(self.valid_component_names) return tuple( - {"expression": eq.evaluate_expression(self, return_type="math_string")} + { + "expression": eq.evaluate_expression( + self, return_type="math_string", references=references + ) + } for eq in equations ) + + @staticmethod + def _dims_to_var_string(da: xr.DataArray) -> str: + if da.shape: + return rf"_\text{{{','.join(str(i).removesuffix('s') for i in da.dims)}}}" + else: + return "" diff --git a/src/calliope/backend/parsing.py b/src/calliope/backend/parsing.py index 53478720..278f037a 100644 --- a/src/calliope/backend/parsing.py +++ b/src/calliope/backend/parsing.py @@ -25,33 +25,45 @@ TRUE_ARRAY = xr.DataArray(True) -class UnparsedEquationDict(TypedDict): - """Unparsed equation checker class.""" +class UnparsedEquation(TypedDict): + """Unparsed equation type hint class.""" where: NotRequired[str] expression: str -class UnparsedConstraintDict(TypedDict): - """Unparsed constraint checker class.""" +class UnparsedConstraint(TypedDict): + """Unparsed constraint type hint class.""" description: NotRequired[str] foreach: NotRequired[list] where: NotRequired[str] - equations: Required[list[UnparsedEquationDict]] - sub_expressions: NotRequired[dict[str, list[UnparsedEquationDict]]] - slices: NotRequired[dict[str, list[UnparsedEquationDict]]] + equations: Required[list[UnparsedEquation]] + sub_expressions: NotRequired[dict[str, list[UnparsedEquation]]] + slices: NotRequired[dict[str, list[UnparsedEquation]]] -class UnparsedExpressionDict(UnparsedConstraintDict): - """Unparsed expression checker class.""" +class UnparsedPiecewiseConstraint(TypedDict): + """Unparsed piecewise constraint type hint class.""" + + description: NotRequired[str] + foreach: NotRequired[list] + where: NotRequired[str] + x_expression: Required[str] + x_values: Required[str] + y_expression: Required[str] + y_values: Required[str] + + +class UnparsedExpression(UnparsedConstraint): + """Unparsed expression type hint class.""" title: NotRequired[str] unit: NotRequired[str] -class UnparsedVariableBoundDict(TypedDict): - """Unparsed variable bounds checker class.""" +class UnparsedVariableBound(TypedDict): + """Unparsed variable bounds type hint class.""" min: str max: str @@ -59,7 +71,7 @@ class UnparsedVariableBoundDict(TypedDict): scale: NotRequired[str] -class UnparsedVariableDict(TypedDict): +class UnparsedVariable(TypedDict): """Unparsed variable checker class.""" title: NotRequired[str] @@ -68,23 +80,24 @@ class UnparsedVariableDict(TypedDict): foreach: list[str] where: str domain: NotRequired[str] - bounds: UnparsedVariableBoundDict + bounds: UnparsedVariableBound -class UnparsedObjectiveDict(TypedDict): +class UnparsedObjective(TypedDict): """Unparsed model objective checker.""" description: NotRequired[str] - equations: Required[list[UnparsedEquationDict]] - sub_expressions: NotRequired[dict[str, list[UnparsedEquationDict]]] + equations: Required[list[UnparsedEquation]] + sub_expressions: NotRequired[dict[str, list[UnparsedEquation]]] sense: str UNPARSED_DICTS = ( - UnparsedConstraintDict - | UnparsedVariableDict - | UnparsedExpressionDict - | UnparsedObjectiveDict + UnparsedConstraint + | UnparsedVariable + | UnparsedExpression + | UnparsedObjective + | UnparsedPiecewiseConstraint ) T = TypeVar("T", bound=UNPARSED_DICTS) @@ -331,7 +344,7 @@ def drop_dims_not_in_foreach(self, where: xr.DataArray) -> xr.DataArray: @overload def evaluate_expression( self, - backend_interface: backend_model.BackendModel, + backend_interface: backend_model.BackendModelGenerator, *, return_type: Literal["array"] = "array", references: set | None = None, @@ -342,7 +355,7 @@ def evaluate_expression( @overload def evaluate_expression( self, - backend_interface: backend_model.BackendModel, + backend_interface: backend_model.BackendModelGenerator, *, return_type: Literal["math_string"], references: set | None = None, @@ -350,7 +363,7 @@ def evaluate_expression( def evaluate_expression( self, - backend_interface: backend_model.BackendModel, + backend_interface: backend_model.BackendModelGenerator, *, return_type: Literal["array", "math_string"] = "array", references: set | None = None, @@ -448,11 +461,18 @@ class ParsedBackendComponent(ParsedBackendEquation): "global_expressions": expression_parser.generate_arithmetic_parser, "objectives": expression_parser.generate_arithmetic_parser, "variables": lambda x: None, + "piecewise_constraints": expression_parser.generate_arithmetic_parser, } def __init__( self, - group: Literal["variables", "global_expressions", "constraints", "objectives"], + group: Literal[ + "variables", + "global_expressions", + "constraints", + "piecewise_constraints", + "objectives", + ], name: str, unparsed_data: T, ) -> None: @@ -538,7 +558,7 @@ def parse_equations( List of parsed equations ready to be evaluated. The length of the list depends on the product of provided equations and sub-expression/slice references. """ - equation_expression_list: list[UnparsedEquationDict] + equation_expression_list: list[UnparsedEquation] equation_expression_list = self._unparsed.get("equations", []) equations = self.generate_expression_list( @@ -639,7 +659,7 @@ def parse_where_string(self, where_string: str = "True") -> pp.ParseResults: def generate_expression_list( self, expression_parser: pp.ParserElement, - expression_list: list[UnparsedEquationDict], + expression_list: list[UnparsedEquation], expression_group: Literal["equations", "sub_expressions", "slices"], id_prefix: str = "", ) -> list[ParsedBackendEquation]: @@ -650,7 +670,7 @@ def generate_expression_list( Args: expression_parser (pp.ParserElement): parser to use. - expression_list (list[UnparsedEquationDict]): list of constraint equations + expression_list (list[UnparsedEquation]): list of constraint equations or sub-expressions with arithmetic expression string and optional where string. expression_group (Literal["equations", "sub_expressions", "slices"]): diff --git a/src/calliope/backend/pyomo_backend_model.py b/src/calliope/backend/pyomo_backend_model.py index c0ca7860..8547f578 100644 --- a/src/calliope/backend/pyomo_backend_model.py +++ b/src/calliope/backend/pyomo_backend_model.py @@ -18,6 +18,11 @@ import pyomo.kernel as pmo # type: ignore import xarray as xr from pyomo.common.tempfiles import TempfileManager # type: ignore +from pyomo.core.kernel.piecewise_library.transforms import ( + PiecewiseLinearFunction, + PiecewiseValidationError, + piecewise_sos2, +) from pyomo.opt import SolverFactory # type: ignore from pyomo.util.model_size import build_model_size_report # type: ignore @@ -28,7 +33,12 @@ T = TypeVar("T") _COMPONENTS_T = Literal[ - "variables", "constraints", "objectives", "parameters", "global_expressions" + "variables", + "constraints", + "piecewise_constraints", + "objectives", + "parameters", + "global_expressions", ] LOGGER = logging.getLogger(__name__) @@ -38,6 +48,7 @@ "variable": "variable", "global_expression": "expression", "constraint": "constraint", + "piecewise_constraint": "block", "objective": "objective", } @@ -58,6 +69,7 @@ def __init__(self, inputs: xr.Dataset, **kwargs) -> None: self._instance.variables = pmo.variable_dict() self._instance.global_expressions = pmo.expression_dict() self._instance.constraints = pmo.constraint_dict() + self._instance.piecewise_constraints = pmo.block_dict() self._instance.objectives = pmo.objective_dict() self._instance.dual = pmo.suffix(direction=pmo.suffix.IMPORT) @@ -69,37 +81,35 @@ def add_parameter( # noqa: D102, override self, parameter_name: str, parameter_values: xr.DataArray, default: Any = np.nan ) -> None: self._raise_error_on_preexistence(parameter_name, "parameters") - - self._create_obj_list(parameter_name, "parameters") - - parameter_da = self._apply_func( - self._to_pyomo_param, - parameter_values.notnull(), - 1, - parameter_values, - name=parameter_name, - ) - - if parameter_da.isnull().all(): + if parameter_values.isnull().all(): self.log( "parameters", parameter_name, "Component not added; no data found in array.", ) - self.delete_component(parameter_name, "parameters") - parameter_da = parameter_da.astype(float) + parameter_da = parameter_values.astype(float) + else: + self._create_obj_list(parameter_name, "parameters") + parameter_da = self._apply_func( + self._to_pyomo_param, + parameter_values.notnull(), + 1, + parameter_values, + name=parameter_name, + ) - parameter_da.attrs["original_dtype"] = parameter_values.dtype attrs = { "title": self._PARAM_TITLES.get(parameter_name, None), "description": self._PARAM_DESCRIPTIONS.get(parameter_name, None), "unit": self._PARAM_UNITS.get(parameter_name, None), "default": default, + "original_dtype": parameter_values.dtype.name, } + self._add_to_dataset(parameter_name, parameter_da, "parameters", attrs) def add_constraint( # noqa: D102, override - self, name: str, constraint_dict: parsing.UnparsedConstraintDict | None = None + self, name: str, constraint_dict: parsing.UnparsedConstraint ) -> None: def _constraint_setter( element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set @@ -122,7 +132,7 @@ def _constraint_setter( self._add_component(name, constraint_dict, _constraint_setter, "constraints") def add_global_expression( # noqa: D102, override - self, name: str, expression_dict: parsing.UnparsedExpressionDict | None = None + self, name: str, expression_dict: parsing.UnparsedExpression ) -> None: def _expression_setter( element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set @@ -139,11 +149,9 @@ def _expression_setter( ) def add_variable( # noqa: D102, override - self, name: str, variable_dict: parsing.UnparsedVariableDict | None = None + self, name: str, variable_dict: parsing.UnparsedVariable ) -> None: domain_dict = {"real": pmo.RealSet, "integer": pmo.IntegerSet} - if variable_dict is None: - variable_dict = self.inputs.attrs["math"]["variables"][name] def _variable_setter(where, references): domain_type = domain_dict[variable_dict.get("domain", "real")] @@ -165,12 +173,10 @@ def _variable_setter(where, references): self._add_component(name, variable_dict, _variable_setter, "variables") def add_objective( # noqa: D102, override - self, name: str, objective_dict: parsing.UnparsedObjectiveDict | None = None + self, name: str, objective_dict: parsing.UnparsedObjective ) -> None: sense_dict = {"minimize": 1, "minimise": 1, "maximize": -1, "maximise": -1} - if objective_dict is None: - objective_dict = self.inputs.attrs["math"]["objectives"][name] sense = sense_dict[objective_dict["sense"]] def _objective_setter( @@ -194,26 +200,17 @@ def _objective_setter( def get_parameter( # noqa: D102, override self, name: str, as_backend_objs: bool = True ) -> xr.DataArray: - parameter = self.parameters.get(name, None) - if parameter is None: - raise KeyError(f"Unknown parameter: {name}") + parameter = self._get_component(name, "parameters") - if as_backend_objs or not isinstance(parameter, xr.DataArray): + if as_backend_objs: return parameter param_as_vals = self._apply_func( self._from_pyomo_param, parameter.notnull(), 1, parameter ) - if parameter.original_dtype.kind == "M": # i.e., np.datetime64 - self.log("parameters", name, "Converting Pyomo object to datetime dtype.") - return xr.apply_ufunc(pd.to_datetime, param_as_vals) - else: - self.log( - "parameters", - name, - f"Converting Pyomo object to {parameter.original_dtype} dtype.", - ) - return param_as_vals.astype(parameter.original_dtype) + orig_dtype = parameter.original_dtype + self.log("parameters", name, f"Converting Pyomo object to {orig_dtype} dtype.") + return param_as_vals.astype(orig_dtype) @overload def get_constraint( # noqa: D102, override @@ -244,9 +241,8 @@ def get_constraint( # noqa: D102, override def get_variable( # noqa: D102, override self, name: str, as_backend_objs: bool = True ) -> xr.DataArray: - variable = self.variables.get(name, None) - if variable is None: - raise KeyError(f"Unknown variable: {name}") + variable = self._get_component(name, "variables") + if as_backend_objs: return variable else: @@ -338,7 +334,12 @@ def __renamer(val, *idx): val.calliope_coords = idx with self._datetime_as_string(self._dataset): - for component_type in ["parameters", "variables", "constraints"]: + for component_type in [ + "parameters", + "variables", + "constraints", + "piecewise_constraints", + ]: for da in self._dataset.filter_by_attrs( coords_in_name=False, **{"obj_type": component_type} ).values(): @@ -372,7 +373,9 @@ def _create_obj_list(self, key: str, component_type: _COMPONENTS_T) -> None: pmo, f"{COMPONENT_TRANSLATOR[singular_component]}_list" )() - def delete_component(self, key: str, component_type: _COMPONENTS_T) -> None: # noqa: D102, override + def delete_component( # noqa: D102, override + self, key: str, component_type: _COMPONENTS_T + ) -> None: component_dict = getattr(self._instance, component_type) if key in component_dict: del component_dict[key] @@ -503,15 +506,33 @@ def unfix_variable( # noqa: D102, override self._apply_func(self._unfix_pyomo_variable, where_da, 1, variable_da) @property - def has_integer_or_binary_variables( # noqa: D102, override - self, - ) -> bool: + def has_integer_or_binary_variables(self) -> bool: # noqa: D102, override model_report = build_model_size_report(self._instance) binaries = model_report["activated"]["binary_variables"] integers = model_report["activated"]["integer_variables"] number_of_binary_and_integer_vars = binaries + integers return number_of_binary_and_integer_vars > 0 + def _to_piecewise_constraint( # noqa: D102, override + self, x_var: Any, y_var: Any, *vals: float, name: str, n_breakpoints: int + ) -> type[ObjPiecewiseConstraint]: + y_vals = pd.Series(vals[n_breakpoints:]).dropna() + x_vals = pd.Series(vals[:n_breakpoints]).dropna() + try: + var = ObjPiecewiseConstraint( + breakpoints=x_vals, + values=y_vals, + input=x_var, + output=y_var, + require_bounded_input_variable=False, + ) + self._instance.piecewise_constraints[name].append(var) + except (PiecewiseValidationError, ValueError) as err: + # We don't want to confuse the user with suggestions of pyomo options they can't access. + err_message = err.args[0].split(" To avoid this error")[0] + raise BackendError(err_message) + return var + def _to_pyomo_param(self, val: Any, *, name: str) -> type[ObjParameter] | float: """Utility function to generate a pyomo parameter for every element of an xarray DataArray. @@ -814,6 +835,26 @@ def getname(self, *args, **kwargs): return self._update_name(pmo.constraint.getname(self, *args, **kwargs)) +class ObjPiecewiseConstraint(piecewise_sos2, CoordObj): + """Pyomo SOS2 piecewise constraint wrapper.""" + + def __init__(self, **kwds): + """Create a Pyomo SOS2 piecesise constraint. + + Created with a `name` property setter (via the `piecewise_sos2.getname` method) + which replaces a list position as a name with a list of strings. + """ + func = PiecewiseLinearFunction( + breakpoints=kwds.pop("breakpoints"), values=kwds.pop("values") + ) + piecewise_sos2.__init__(self, func, **kwds) + CoordObj.__init__(self) + + def getname(self, *args, **kwargs): + """Get piecewise constraint name.""" + return self._update_name(piecewise_sos2.getname(self, *args, **kwargs)) + + class PyomoShadowPrices(backend_model.ShadowPrices): """Pyomo shadow price functionality.""" diff --git a/src/calliope/backend/where_parser.py b/src/calliope/backend/where_parser.py index 110132e4..f434a9bf 100644 --- a/src/calliope/backend/where_parser.py +++ b/src/calliope/backend/where_parser.py @@ -153,7 +153,7 @@ def __repr__(self): """Programming / official string representation.""" return f"DATA_VAR:{self.data_var}" - def _preprocess(self) -> tuple[xr.Dataset, str]: + def _preprocess(self) -> str: """Get data variable from the optimisation problem dataset. Raises: @@ -179,13 +179,7 @@ def _preprocess(self) -> tuple[xr.Dataset, str]: "These arrays cannot be used for comparison with expected values. " f"Received `{self.data_var}`." ) - - if data_var_type == "parameters": - source_dataset = self.eval_attrs["input_data"] - else: - source_dataset = backend_interface._dataset - - return source_dataset, data_var_type + return data_var_type def _data_var_exists( self, source_dataset: xr.Dataset, data_var_type: str @@ -212,24 +206,26 @@ def _data_var_with_default(self, source_dataset: xr.Dataset) -> xr.DataArray: return var def as_math_string(self) -> str: # noqa: D102, override - # TODO: add dims from a YAML schema of params that includes default dims - source_dataset, data_var_type = self._preprocess() - if data_var_type == "parameters": - data_var_string = rf"\textit{{{self.data_var}}}" - else: - data_var_string = rf"\textbf{{{self.data_var}}}" + self._preprocess() - var = source_dataset.get(self.data_var, None) - if var is not None and var.shape: - data_var_string += ( - rf"_\text{{{','.join(str(i).removesuffix('s') for i in var.dims)}}}" - ) + var = self.eval_attrs["backend_interface"]._dataset.get( + self.data_var, xr.DataArray() + ) + + try: + data_var_string = var.attrs["math_repr"] + except (AttributeError, KeyError): + data_var_string = rf"\text{{{self.data_var}}}" if self.eval_attrs.get("apply_where", True): data_var_string = rf"\exists ({data_var_string})" return data_var_string def as_array(self) -> xr.DataArray: # noqa: D102, override - source_dataset, data_var_type = self._preprocess() + data_var_type = self._preprocess() + if data_var_type == "parameters": + source_dataset = self.eval_attrs["input_data"] + else: + source_dataset = self.eval_attrs["backend_interface"]._dataset if self.eval_attrs.get("apply_where", True): return self._data_var_exists(source_dataset, data_var_type) diff --git a/src/calliope/config/math_schema.yaml b/src/calliope/config/math_schema.yaml index a75b3602..38310d81 100644 --- a/src/calliope/config/math_schema.yaml +++ b/src/calliope/config/math_schema.yaml @@ -86,6 +86,30 @@ properties: type: string description: Index slice expression, such as a list of set items or a call to a helper function. + piecewise_constraints: + type: object + description: All _piecewise_ constraints to apply to the optimisation problem. + additionalProperties: false + patternProperties: + '[^\d^_\W][\w\d]+': + type: object + description: A named piecewise constraint, linking an `x`-axis decision variable with a `y`-axis decision variable with values at specified breakpoints. + additionalProperties: false + required: [x_expression, x_values, y_expression, y_values] + properties: + description: *description + active: *active + foreach: *foreach + where: *top_level_where + x_expression: &piecewisevar + type: string + description: Variable name whose values are assigned at each breakpoint. + x_values: &piecewisevals + type: string + description: Parameter name containing data, indexed over the `breakpoints` dimension. + y_expression: *piecewisevar + y_values: *piecewisevals + global_expressions: type: object description: >- diff --git a/src/calliope/math/base.yaml b/src/calliope/math/base.yaml index 36ba1daf..1d9d6fc4 100644 --- a/src/calliope/math/base.yaml +++ b/src/calliope/math/base.yaml @@ -930,4 +930,6 @@ global_expressions: - where: "NOT cost_var" expression: "0" active: true # optional; defaults to true. - # --8<-- [end:expression] \ No newline at end of file + # --8<-- [end:expression] + +piecewise_constraints: {} \ No newline at end of file diff --git a/tests/common/lp_files/sos2_piecewise_cost_investment.lp b/tests/common/lp_files/sos2_piecewise_cost_investment.lp new file mode 100644 index 00000000..7a7716ae --- /dev/null +++ b/tests/common/lp_files/sos2_piecewise_cost_investment.lp @@ -0,0 +1,46 @@ +\* Source Pyomo model name=None *\ + +min +objectives(dummy_obj)(0): ++2.0 ONE_VAR_CONSTANT + +s.t. + +c_e_piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_c(0)_: ++1.0 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(1) ++2.0 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(2) ++10.0 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(3) +-1 variables(flow_cap)(a__test_supply_elec__electricity) += 0 + +c_e_piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_c(1)_: ++5.0 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(1) ++8.0 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(2) ++20.0 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(3) +-1 variables(piecewise_cost_investment)(a__test_supply_elec__electricity__monetary) += 0 + +c_e_piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_c(2)_: ++1 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(1) ++1 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(2) ++1 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(3) ++1 piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(0) += 1 + +bounds + 1 <= ONE_VAR_CONSTANT <= 1 + 0 <= piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(1) <= +inf + 0 <= piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(2) <= +inf + 0 <= piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(3) <= +inf + 0 <= variables(flow_cap)(a__test_supply_elec__electricity) <= 10.0 + 0 <= variables(piecewise_cost_investment)(a__test_supply_elec__electricity__monetary) <= +inf + 0 <= piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(0) <= +inf +SOS + +piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_s: S2:: + piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(0):1 + piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(1):2 + piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(2):3 + piecewise_constraints(sos2_piecewise_costs)(a__test_supply_elec__electricity__monetary)_v(3):4 + +end diff --git a/tests/common/util.py b/tests/common/util.py index 90d3803a..290081f7 100644 --- a/tests/common/util.py +++ b/tests/common/util.py @@ -107,7 +107,8 @@ def build_lp( getattr(backend_instance, f"add_{component}")(name, dict_) elif isinstance(component_math, list): for name in component_math: - getattr(backend_instance, f"add_{component}")(name) + dict_ = model.math[component_group][name] + getattr(backend_instance, f"add_{component}")(name, dict_) # MUST have an objective for a valid LP file if math is None or "objectives" not in math.keys(): diff --git a/tests/test_backend_general.py b/tests/test_backend_general.py index b26a7c4c..7bd73d11 100644 --- a/tests/test_backend_general.py +++ b/tests/test_backend_general.py @@ -174,7 +174,7 @@ def test_get_parameter(self, parameter): assert parameter.attrs == { "obj_type": "parameters", "is_result": 0, - "original_dtype": np.dtype("float64"), + "original_dtype": "float64", "references": {"flow_in_inc_eff"}, "coords_in_name": False, "default": 1.0, @@ -254,7 +254,24 @@ def test_get_global_expression_as_vals_no_solve(self, built_model_cls_longnames) expr = built_model_cls_longnames.backend.get_global_expression( "cost", as_backend_objs=False, eval_body=True ) - assert expr.to_series().dropna().apply(lambda x: isinstance(x, str)).all() + assert ( + expr.where(expr != "nan") + .to_series() + .dropna() + .apply(lambda x: isinstance(x, str)) + .all() + ) + + def test_timeseries_dtype(self, built_model_cls_longnames): + """Getting verbose strings leads to the timeseries being stringified then converted back to datetime.""" + expr = built_model_cls_longnames.backend.get_global_expression( + "flow_out_inc_eff", as_backend_objs=False, eval_body=True + ) + assert ( + expr.where(expr != "nan").to_series().dropna().str.contains("2005-").all() + ) + assert built_model_cls_longnames.backend._dataset.timesteps.dtype.kind == "M" + assert built_model_cls_longnames.backend.inputs.timesteps.dtype.kind == "M" def test_get_constraint_attrs(self, constraint): """Check a constraint has all expected attributes.""" @@ -675,3 +692,109 @@ def test_has_integer_or_binary_variables_lp(self, solved_model_cls): def test_has_integer_or_binary_variables_milp(self, solved_model_milp_cls): """MILP models have integer / binary variables.""" assert solved_model_milp_cls.backend.has_integer_or_binary_variables + + +class TestPiecewiseConstraints: + def gen_params(self, data, index=[0, 1, 2], dim="breakpoints"): + return { + "parameters": { + "piecewise_x": {"data": data, "index": index, "dims": dim}, + "piecewise_y": { + "data": [0, 1, 5], + "index": [0, 1, 2], + "dims": "breakpoints", + }, + } + } + + @pytest.fixture(scope="class") + def working_math(self): + return { + "foreach": ["nodes", "techs", "carriers"], + "where": "[test_supply_elec] in techs AND piecewise_x AND piecewise_y", + "x_values": "piecewise_x", + "x_expression": "flow_cap", + "y_values": "piecewise_y", + "y_expression": "source_cap", + "description": "FOO", + } + + @pytest.fixture(scope="class") + def working_params(self): + return self.gen_params([0, 5, 10]) + + @pytest.fixture(scope="class") + def length_mismatch_params(self): + return self.gen_params([0, 10], [0, 1]) + + @pytest.fixture(scope="class") + def not_reaching_var_bound_with_breakpoint_params(self): + return self.gen_params([0, 5, 8]) + + @pytest.fixture(scope="class") + def missing_breakpoint_dims(self): + return self.gen_params([0, 5, 10], dim="foobar") + + @pytest.fixture(scope="class") + def working_model(self, backend, working_params, working_math): + m = build_model(working_params, "simple_supply,two_hours,investment_costs") + m.build(backend=backend) + m.backend.add_piecewise_constraint("foo", working_math) + return m + + @pytest.fixture(scope="class") + def piecewise_constraint(self, working_model): + return working_model.backend.get_piecewise_constraint("foo") + + def test_piecewise_attrs(self, piecewise_constraint): + """Check a piecewise constraint has all expected attributes.""" + expected_keys = set( + ["obj_type", "references", "description", "yaml_snippet", "coords_in_name"] + ) + assert not expected_keys.symmetric_difference(piecewise_constraint.attrs.keys()) + + def test_piecewise_obj_type(self, piecewise_constraint): + """Check a piecewise constraint has expected object type.""" + assert piecewise_constraint.attrs["obj_type"] == "piecewise_constraints" + + def test_piecewise_refs(self, piecewise_constraint): + """Check a piecewise constraint has expected refs to other math components (zero for piecewise constraints).""" + assert not piecewise_constraint.attrs["references"] + + def test_piecewise_obj_coords_in_name(self, piecewise_constraint): + """Check a piecewise constraint does not have verbose strings activated.""" + assert piecewise_constraint.attrs["coords_in_name"] is False + + @pytest.mark.parametrize( + "var", ["flow_cap", "source_cap", "piecewise_x", "piecewise_y"] + ) + def test_piecewise_upstream_refs(self, working_model, var): + """Expected tracking of piecewise constraint in component reference chains.""" + assert "foo" in working_model.backend._dataset[var].attrs["references"] + + def test_fails_on_breakpoints_in_foreach(self, working_model, working_math): + """Expected error when defining `breakpoints` in foreach.""" + failing_math = {"foreach": ["nodes", "techs", "carriers", "breakpoints"]} + with pytest.raises(calliope.exceptions.BackendError) as excinfo: + working_model.backend.add_piecewise_constraint( + "bar", {**working_math, **failing_math} + ) + assert check_error_or_warning( + excinfo, + "(piecewise_constraints, bar) | `breakpoints` dimension should not be in `foreach`", + ) + + def test_fails_on_no_breakpoints_in_params( + self, missing_breakpoint_dims, working_math, backend + ): + """Expected error when parameter defining breakpoints isn't indexed over `breakpoints`.""" + m = build_model( + missing_breakpoint_dims, "simple_supply,two_hours,investment_costs" + ) + m.build(backend=backend) + with pytest.raises(calliope.exceptions.BackendError) as excinfo: + m.backend.add_piecewise_constraint("bar", working_math) + assert check_error_or_warning( + excinfo, + "(piecewise_constraints, bar) | `x_values` must be indexed over the `breakpoints` dimension", + ) diff --git a/tests/test_backend_gurobi.py b/tests/test_backend_gurobi.py index 4ac7a8b0..3c0fef86 100755 --- a/tests/test_backend_gurobi.py +++ b/tests/test_backend_gurobi.py @@ -302,3 +302,93 @@ def test_get_shadow_price_empty_milp(self, supply_milp): supply_milp.solve() shadow_prices = supply_milp.backend.shadow_prices.get("system_balance") assert shadow_prices.isnull().all() + + +class TestPiecewiseConstraints: + def gen_params(self, data, index=[0, 1, 2], dim="breakpoints"): + return { + "parameters": { + "piecewise_x": {"data": data, "index": index, "dims": dim}, + "piecewise_y": { + "data": [0, 1, 5], + "index": [0, 1, 2], + "dims": "breakpoints", + }, + } + } + + @pytest.fixture(scope="class") + def working_math(self): + return { + "foreach": ["nodes", "techs", "carriers"], + "where": "[test_supply_elec] in techs AND piecewise_x AND piecewise_y", + "x_values": "piecewise_x", + "x_expression": "flow_cap", + "y_values": "piecewise_y", + "y_expression": "source_cap", + "description": "FOO", + } + + @pytest.fixture(scope="class") + def failing_math(self, working_math): + return {**working_math, **{"y_expression": "sum(flow_in, over=timesteps)"}} + + @pytest.fixture(scope="class") + def working_params(self): + return self.gen_params([0, 5, 10]) + + @pytest.fixture(scope="class") + def length_mismatch_params(self): + return self.gen_params([0, 10], [0, 1]) + + @pytest.fixture(scope="class") + def working_model(self, working_params, working_math): + m = build_model(working_params, "simple_supply,two_hours,investment_costs") + m.build(backend="gurobi") + m.backend.add_piecewise_constraint("foo", working_math) + return m + + def test_piecewise_type(self, working_model): + """All piecewise elements are the correct Gurobi type.""" + constr = working_model.backend.get_piecewise_constraint("foo") + assert ( + constr.to_series() + .dropna() + .apply(lambda x: isinstance(x, gurobipy.GenConstr)) + .all() + ) + + def test_piecewise_verbose(self, working_model): + """All piecewise elements have the full set of dimensions when verbose.""" + working_model.backend.verbose_strings() + constr = working_model.backend.get_piecewise_constraint("foo") + dims = {"nodes": "a", "techs": "test_supply_elec", "carriers": "electricity"} + constraint_item = constr.sel(dims).item() + assert ( + constraint_item.GenConstrName + == f"foo[{', '.join(dims[i] for i in constr.dims)}]" + ) + + def test_fails_on_length_mismatch(self, length_mismatch_params, working_math): + """Expected error when number of breakpoints on X and Y don't match.""" + m = build_model( + length_mismatch_params, "simple_supply,two_hours,investment_costs" + ) + m.build(backend="gurobi") + with pytest.raises(exceptions.BackendError) as excinfo: + m.backend.add_piecewise_constraint("foo", working_math) + assert check_error_or_warning( + excinfo, + "Errors in generating piecewise constraint: Arguments xpts and ypts must have the same length", + ) + + def test_expressions_not_allowed(self, working_params, failing_math): + """Expected error when using an expression instead of a decision variable (gurobi-specific error).""" + m = build_model(working_params, "simple_supply,two_hours,investment_costs") + m.build(backend="gurobi") + with pytest.raises(exceptions.BackendError) as excinfo: + m.backend.add_piecewise_constraint("foo", failing_math) + assert check_error_or_warning( + excinfo, + "Gurobi backend can only build piecewise constraints using decision variables.", + ) diff --git a/tests/test_backend_latex_backend.py b/tests/test_backend_latex_backend.py index 69713445..f1bac0a0 100644 --- a/tests/test_backend_latex_backend.py +++ b/tests/test_backend_latex_backend.py @@ -68,6 +68,9 @@ def test_invalid_format(self, build_all, tmpdir_factory, filepath, format): class TestLatexBackendModel: + def test_inputs(self, dummy_latex_backend_model, dummy_model_data): + assert dummy_latex_backend_model.inputs.equals(dummy_model_data) + @pytest.mark.parametrize( "backend_obj", ["valid_latex_backend", "dummy_latex_backend_model"] ) @@ -97,6 +100,10 @@ def test_add_variable(self, request, dummy_model_data, backend_obj): ) assert "var" in latex_backend_model.valid_component_names assert "math_string" in latex_backend_model.variables["var"].attrs + assert ( + latex_backend_model.variables["var"].attrs["math_repr"] + == r"\textbf{var}_\text{node,tech}" + ) def test_add_variable_not_valid(self, valid_latex_backend): valid_latex_backend.add_variable( @@ -249,6 +256,79 @@ def test_add_objective(self, dummy_latex_backend_model): assert "obj" not in dummy_latex_backend_model.valid_component_names assert len(dummy_latex_backend_model.objectives.data_vars) == 1 + def test_add_piecewise_constraint(self, dummy_latex_backend_model): + dummy_latex_backend_model.add_parameter( + "piecewise_x", + xr.DataArray(data=[0, 5, 10], coords={"breakpoints": [0, 1, 2]}), + ) + dummy_latex_backend_model.add_parameter( + "piecewise_y", + xr.DataArray(data=[0, 1, 5], coords={"breakpoints": [0, 1, 2]}), + ) + for param in ["piecewise_x", "piecewise_y"]: + dummy_latex_backend_model.inputs[param] = ( + dummy_latex_backend_model._dataset[param] + ) + dummy_latex_backend_model.add_piecewise_constraint( + "p_constr", + { + "foreach": ["nodes", "techs"], + "where": "piecewise_x AND piecewise_y", + "x_values": "piecewise_x", + "x_expression": "multi_dim_var + 1", + "y_values": "piecewise_y", + "y_expression": "no_dim_var", + "description": "FOO", + }, + ) + math_string = dummy_latex_backend_model.piecewise_constraints["p_constr"].attrs[ + "math_string" + ] + assert ( + r"\text{ breakpoint }\negthickspace \in \negthickspace\text{ breakpoints }" + in math_string + ) + assert ( + r"\text{if } \textbf{multi_dim_var}_\text{node,tech} + 1\mathord{=}\textit{piecewise_x}_\text{breakpoint}" + in math_string + ) + assert ( + r"\textbf{no_dim_var}\mathord{=}\textit{piecewise_y}_\text{breakpoint}" + in math_string + ) + + def test_add_piecewise_constraint_no_foreach(self, dummy_latex_backend_model): + dummy_latex_backend_model.add_parameter( + "piecewise_x", + xr.DataArray(data=[0, 5, 10], coords={"breakpoints": [0, 1, 2]}), + ) + dummy_latex_backend_model.add_parameter( + "piecewise_y", + xr.DataArray(data=[0, 1, 5], coords={"breakpoints": [0, 1, 2]}), + ) + for param in ["piecewise_x", "piecewise_y"]: + dummy_latex_backend_model.inputs[param] = ( + dummy_latex_backend_model._dataset[param] + ) + dummy_latex_backend_model.add_piecewise_constraint( + "p_constr_no_foreach", + { + "where": "piecewise_x AND piecewise_y", + "x_values": "piecewise_x", + "x_expression": "sum(multi_dim_var, over=[nodes, techs])", + "y_values": "piecewise_y", + "y_expression": "no_dim_var", + "description": "BAR", + }, + ) + math_string = dummy_latex_backend_model.piecewise_constraints[ + "p_constr_no_foreach" + ].attrs["math_string"] + assert ( + r"\text{ breakpoint }\negthickspace \in \negthickspace\text{ breakpoints }" + in math_string + ) + def test_create_obj_list(self, dummy_latex_backend_model): assert dummy_latex_backend_model._create_obj_list("var", "variables") is None @@ -555,10 +635,12 @@ def test_render(self, dummy_latex_backend_model, instring, kwargs, expected): def test_get_variable_bounds_string(self, dummy_latex_backend_model): bounds = {"min": 1, "max": 2e6} + refs = set() lb, ub = dummy_latex_backend_model._get_variable_bounds_string( - "multi_dim_var", bounds + "multi_dim_var", bounds, refs ) assert lb == {"expression": r"1 \leq \textbf{multi_dim_var}_\text{node,tech}"} assert ub == { "expression": r"\textbf{multi_dim_var}_\text{node,tech} \leq 2\mathord{\times}10^{+06}" } + assert refs == {"multi_dim_var"} diff --git a/tests/test_backend_pyomo.py b/tests/test_backend_pyomo.py index 08780750..61bf37f6 100755 --- a/tests/test_backend_pyomo.py +++ b/tests/test_backend_pyomo.py @@ -12,6 +12,7 @@ import xarray as xr from calliope.attrdict import AttrDict from calliope.backend.pyomo_backend_model import PyomoBackendModel +from pyomo.core.kernel.piecewise_library.transforms import piecewise_sos2 from .common.util import build_test_model as build_model from .common.util import check_error_or_warning, check_variable_exists @@ -1692,6 +1693,22 @@ def test_new_build_get_variable_as_vals(self, simple_supply): .any() ) + def test_new_build_get_component_exists(self, simple_supply): + param = simple_supply.backend._get_component("flow_in_eff", "parameters") + assert isinstance(param, xr.DataArray) + + def test_new_build_get_component_does_not_exist(self, simple_supply): + with pytest.raises(KeyError) as excinfo: + simple_supply.backend._get_component("does_not_exist", "parameters") + assert check_error_or_warning(excinfo, "Unknown parameter: does_not_exist") + + def test_new_build_get_component_wrong_group(self, simple_supply): + with pytest.raises(KeyError) as excinfo: + simple_supply.backend._get_component("flow_in_eff", "piecewise_constraints") + assert check_error_or_warning( + excinfo, "Unknown piecewise constraint: flow_in_eff" + ) + def test_new_build_get_parameter(self, simple_supply): """Check a parameter has the correct data type and has all expected attributes.""" param = simple_supply.backend.get_parameter("cost_flow_cap") @@ -2018,6 +2035,105 @@ def test_verbose_strings_no_len(self, simple_supply_longnames): assert obj.coords_in_name +class TestPiecewiseConstraints: + def gen_params(self, data, index=[0, 1, 2], dim="breakpoints"): + return { + "parameters": { + "piecewise_x": {"data": data, "index": index, "dims": dim}, + "piecewise_y": { + "data": [0, 1, 5], + "index": [0, 1, 2], + "dims": "breakpoints", + }, + } + } + + @pytest.fixture(scope="class") + def working_math(self): + return { + "foreach": ["nodes", "techs", "carriers"], + "where": "[test_supply_elec] in techs AND piecewise_x AND piecewise_y", + "x_values": "piecewise_x", + "x_expression": "flow_cap", + "y_values": "piecewise_y", + "y_expression": "sum(flow_in, over=timesteps)", + "description": "FOO", + } + + @pytest.fixture(scope="class") + def working_params(self): + return self.gen_params([0, 5, 10]) + + @pytest.fixture(scope="class") + def length_mismatch_params(self): + return self.gen_params([0, 10], [0, 1]) + + @pytest.fixture(scope="class") + def not_reaching_var_bound_with_breakpoint_params(self): + return self.gen_params([0, 5, 8]) + + @pytest.fixture(scope="class") + def working_model(self, working_params, working_math): + m = build_model(working_params, "simple_supply,two_hours,investment_costs") + m.build() + m.backend.add_piecewise_constraint("foo", working_math) + return m + + def test_piecewise_type(self, working_model): + """All piecewise elements are the correct Pyomo type.""" + constr = working_model.backend.get_piecewise_constraint("foo") + assert ( + constr.to_series() + .dropna() + .apply(lambda x: isinstance(x, piecewise_sos2)) + .all() + ) + + def test_piecewise_verbose(self, working_model): + """All piecewise elements have the full set of dimensions when verbose.""" + working_model.backend.verbose_strings() + constr = working_model.backend.get_piecewise_constraint("foo") + dims = {"nodes": "a", "techs": "test_supply_elec", "carriers": "electricity"} + constraint_item = constr.sel(dims).item() + assert ( + str(constraint_item) + == f"piecewise_constraints[foo][{', '.join(dims[i] for i in constr.dims)}]" + ) + + def test_fails_on_length_mismatch(self, length_mismatch_params, working_math): + """Expected error when number of breakpoints on X and Y don't match.""" + m = build_model( + length_mismatch_params, "simple_supply,two_hours,investment_costs" + ) + m.build() + with pytest.raises(exceptions.BackendError) as excinfo: + m.backend.add_piecewise_constraint("foo", working_math) + assert check_error_or_warning( + excinfo, + "The number of breakpoints (2) differs from the number of function values (3)", + ) + + def test_fails_on_not_reaching_bounds( + self, not_reaching_var_bound_with_breakpoint_params, working_math + ): + """Expected error when breakpoints exceed upper bound of the variable (pyomo-specific error).""" + m = build_model( + not_reaching_var_bound_with_breakpoint_params, + "simple_supply,two_hours,investment_costs", + ) + m.build() + with pytest.raises(exceptions.BackendError) as excinfo: + m.backend.add_piecewise_constraint("foo", working_math) + assert check_error_or_warning( + excinfo, + [ + "(piecewise_constraints, foo) | Errors in generating piecewise constraint: Piecewise function domain does not include the upper bound", + "ub = 10.0 > 8.0.", + ], + ) + assert not check_error_or_warning(excinfo, "To avoid this error") + + class TestShadowPrices: @pytest.fixture() def simple_supply(self): diff --git a/tests/test_backend_where_parser.py b/tests/test_backend_where_parser.py index 1e7ec946..a64caead 100644 --- a/tests/test_backend_where_parser.py +++ b/tests/test_backend_where_parser.py @@ -564,15 +564,16 @@ def test_where_malformed(self, where, instring): class TestAsMathString: @pytest.fixture() - def latex_eval_kwargs(self, eval_kwargs): + def latex_eval_kwargs(self, eval_kwargs, dummy_latex_backend_model): eval_kwargs["return_type"] = "math_string" + eval_kwargs["backend_interface"] = dummy_latex_backend_model return eval_kwargs @pytest.mark.parametrize( ("parser", "instring", "expected"), [ ("data_var", "with_inf", r"\exists (\textit{with_inf}_\text{node,tech})"), - ("data_var", "foo", r"\exists (\textit{foo})"), + ("data_var", "foo", r"\exists (\text{foo})"), ("data_var", "no_dims", r"\exists (\textit{no_dims})"), ("config_option", "config.foo", r"\text{config.foo}"), ("bool_operand", "True", "true"), diff --git a/tests/test_math.py b/tests/test_math.py index a37b07b1..3a5f263e 100644 --- a/tests/test_math.py +++ b/tests/test_math.py @@ -514,6 +514,41 @@ def test_piecewise(self, build_and_compare): ) +class TestSOS2PiecewiseCosts(CustomMathExamples): + YAML_FILEPATH = "sos2_piecewise_linear_costs.yaml" + + def test_piecewise(self, build_and_compare): + overrides = { + "techs.test_supply_elec.lifetime": 10, + "techs.test_supply_elec.cost_interest_rate": { + "data": 0.1, + "index": "monetary", + "dims": "costs", + }, + "techs.test_supply_elec": { + "piecewise_cost_investment_x": { + "data": [0, 1, 2, 10], + "index": [0, 1, 2, 3], + "dims": "breakpoints", + }, + "piecewise_cost_investment_y": { + "data": [0, 5, 8, 20], + "index": [0, 1, 2, 3], + "dims": "breakpoints", + }, + }, + } + build_and_compare( + "sos2_piecewise_cost_investment", + "supply_purchase,two_hours", + overrides, + components={ + "piecewise_constraints": ["sos2_piecewise_costs"], + "variables": ["piecewise_cost_investment"], + }, + ) + + class TestPiecewiseEfficiency(CustomMathExamples): YAML_FILEPATH = "piecewise_linear_efficiency.yaml"