Feature: piecewise constraints

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`.

docs/examples/piecewise_constraints.py

@@ -0,0 +1,270 @@
+# ---
+# 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:
+#
+# <div class="admonition note">
+#     <p class="admonition-title">Note</p>
+#     <p>
+#         We must index our piecewise data over "breakpoints".
+#     </p>
+# </div>
+#
+
+# %%
+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",
+        },
+    }
+}
+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.
+#

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.

docs/user_defined_math/components.md

@@ -87,6 +87,50 @@ 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.
+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:

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

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: