diff --git a/CHANGELOG.md b/CHANGELOG.md
index 92163728..867ed39b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,8 @@
 
 ### User-facing changes
 
+|new| Math has been removed from `model.math`, and can now be accessed via `model.math.data` (#639).
+
 |new| (non-NaN) Default values and data types for parameters appear in math documentation (if they appear in the model definition schema) (#677).
 
 |changed| `data_sources` -> `data_tables` and `data_sources.source` -> `data_tables.data`.
@@ -63,6 +65,12 @@ Parameter titles from the model definition schema will also propagate to the mod
 
 ### Internal changes
 
+|changed| `model._model_def_dict` has been removed.
+
+|new| `CalliopeMath` is a new helper class to handle math additions, including separate methods for pre-defined math, user-defined math and validation checks.
+
+|changed| `MathDocumentation` has been extracted from `Model`/`LatexBackend`, and now is a postprocessing module which can take models as input.
+
 |new| `gurobipy` is a development dependency that will be added as an optional dependency to the conda-forge calliope feedstock recipe.
 
 |changed| Added any new math dicts defined with `calliope.Model.backend.add_[...](...)` to the backend math dict registry stored in `calliope.Model.backend.inputs.attrs["math"]`.
diff --git a/docs/examples/calliope_model_object.py b/docs/examples/calliope_model_object.py
index 9236e697..a4930f43 100644
--- a/docs/examples/calliope_model_object.py
+++ b/docs/examples/calliope_model_object.py
@@ -36,33 +36,10 @@
 # Get information on the model
 print(m.info())
 
-# %% [markdown]
-# ## Model definition dictionary
-#
-# `m._model_def_dict` is a python dictionary that holds all the data from the model definition YAML files, restructured into one dictionary.
-#
-# The underscore before the method indicates that it defaults to being hidden (i.e. you wouldn't see it by trying a tab auto-complete and it isn't documented)
-
-# %%
-m._model_def_dict.keys()
-
-# %% [markdown]
-# `techs` hold only the information about a technology that is specific to that node
-
-# %%
-m._model_def_dict["techs"]["pv"]
-
-# %% [markdown]
-# `nodes` hold only the information about a technology that is specific to that node
-
-# %%
-m._model_def_dict["nodes"]["X2"]["techs"]["pv"]
-
 # %% [markdown]
 # ## Model data
 #
-# `m._model_data` is an xarray Dataset.
-# Like `_model_def_dict` it is a hidden prperty of the Model as you are expected to access the data via the public property `inputs`
+# `m._model_data` is an xarray Dataset, a hidden property of the Model as you are expected to access the data via the public property `inputs`
 
 # %%
 m.inputs
diff --git a/docs/examples/piecewise_constraints.py b/docs/examples/piecewise_constraints.py
index 8d0b66db..064be53d 100644
--- a/docs/examples/piecewise_constraints.py
+++ b/docs/examples/piecewise_constraints.py
@@ -68,22 +68,20 @@
 #
 
 # %%
-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",
-        },
-    }
-}
+new_params = f"""
+  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)
+new_params_as_dict = calliope.AttrDict.from_yaml_string(new_params)
+m = calliope.examples.national_scale(override_dict=new_params_as_dict)
 
 # %%
 m.inputs.capacity_steps
@@ -95,55 +93,48 @@
 # ## 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`.
+# In this example, we need:
+# 1. a new decision variable for investment costs that can take on the value defined by the curve at a given value of `flow_cap`;
+# 1. to link that decision variable to our total cost calculation; and
+# 1. to define the piecewise constraint.
 
 # %%
-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:
+new_math = """
+  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
+  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"
+  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]
 # # Building and checking the optimisation problem
 #
-# With our piecewise constraint defined, we can build our optimisation problem
+# With our piecewise constraint defined, we can build our optimisation problem and inject this new math.
 
 # %%
-m.build()
+new_math_as_dict = calliope.AttrDict.from_yaml_string(new_math)
+m.build(add_math_dict=new_math_as_dict)
 
 # %% [markdown]
 # And we can see that our piecewise constraint exists in the built optimisation problem "backend"
@@ -190,65 +181,6 @@
 )
 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
 #
diff --git a/docs/hooks/dummy_model/model.yaml b/docs/hooks/dummy_model/model.yaml
index 81f4fdd2..e9c1976f 100644
--- a/docs/hooks/dummy_model/model.yaml
+++ b/docs/hooks/dummy_model/model.yaml
@@ -2,8 +2,9 @@ overrides:
   storage_inter_cluster:
     config.init:
       name: inter-cluster storage
-      add_math: ["storage_inter_cluster"]
       time_cluster: cluster_days.csv
+    config.build:
+      add_math: ["storage_inter_cluster"]
 
 config.init.name: base
 
diff --git a/docs/hooks/generate_math_docs.py b/docs/hooks/generate_math_docs.py
index 4849da88..b513ebd2 100644
--- a/docs/hooks/generate_math_docs.py
+++ b/docs/hooks/generate_math_docs.py
@@ -11,6 +11,7 @@
 from mkdocs.structure.files import File
 
 import calliope
+from calliope.postprocess.math_documentation import MathDocumentation
 
 logger = logging.getLogger("mkdocs")
 
@@ -43,31 +44,33 @@ def on_files(files: list, config: dict, **kwargs):
     """Process documentation for pre-defined calliope math files."""
     model_config = calliope.AttrDict.from_yaml(MODEL_PATH)
 
-    base_model = generate_base_math_model()
+    base_documentation = generate_base_math_documentation()
     write_file(
-        "base.yaml",
+        "plan.yaml",
         textwrap.dedent(
             """
         Complete base mathematical formulation for a Calliope model.
         This math is _always_ applied but can be overridden with pre-defined additional math or [your own math][adding-your-own-math-to-a-model].
         """
         ),
-        base_model,
+        base_documentation,
         files,
         config,
     )
 
     for override in model_config["overrides"].keys():
-        custom_model = generate_custom_math_model(base_model, override)
+        custom_documentation = generate_custom_math_documentation(
+            base_documentation, override
+        )
         write_file(
             f"{override}.yaml",
             textwrap.dedent(
                 f"""
-            Pre-defined additional math to apply {custom_model.inputs.attrs['name']} math on top of the [base mathematical formulation][base-math].
+            Pre-defined additional math to apply {custom_documentation.name} math on top of the [base mathematical formulation][base-math].
             This math is _only_ applied if referenced in the `config.init.add_math` list as `{override}`.
             """
             ),
-            custom_model,
+            custom_documentation,
             files,
             config,
         )
@@ -78,7 +81,7 @@ def on_files(files: list, config: dict, **kwargs):
 def write_file(
     filename: str,
     description: str,
-    model: calliope.Model,
+    math_documentation: MathDocumentation,
     files: list[File],
     config: dict,
 ) -> None:
@@ -87,12 +90,10 @@ def write_file(
     Args:
         filename (str): name of produced `.md` file.
         description (str): first paragraph after title.
-        model (calliope.Model): calliope model with the given math.
+        math_documentation (MathDocumentation): calliope math documentation.
         files (list[File]): math files to parse.
         config (dict): documentation configuration.
     """
-    title = model.inputs.attrs["name"] + " math"
-
     output_file = (Path("math") / filename).with_suffix(".md")
     output_full_filepath = Path(TEMPDIR.name) / output_file
     output_full_filepath.parent.mkdir(exist_ok=True, parents=True)
@@ -123,7 +124,8 @@ def write_file(
 
     nav_reference["Pre-defined math"].append(output_file.as_posix())
 
-    math_doc = model.math_documentation.write(format="md", mkdocs_features=True)
+    title = math_documentation.name
+    math_doc = math_documentation.write(format="md", mkdocs_features=True)
     file_to_download = Path("..") / filename
     output_full_filepath.write_text(
         PREPEND_SNIPPET.format(
@@ -136,65 +138,67 @@ def write_file(
     )
 
 
-def generate_base_math_model() -> calliope.Model:
-    """Generate model with documentation for the base math.
-
-    Args:
-        model_config (dict): Calliope model config.
+def generate_base_math_documentation() -> MathDocumentation:
+    """Generate model documentation for the base math.
 
     Returns:
-        calliope.Model: Base math model to use in generating math docs.
+        MathDocumentation: model math documentation with latex backend.
     """
     model = calliope.Model(model_definition=MODEL_PATH)
-    model.math_documentation.build()
-    return model
+    model.build()
+    return MathDocumentation(model)
 
 
-def generate_custom_math_model(
-    base_model: calliope.Model, override: str
-) -> calliope.Model:
-    """Generate model with documentation for a pre-defined math file.
+def generate_custom_math_documentation(
+    base_documentation: MathDocumentation, override: str
+) -> MathDocumentation:
+    """Generate model documentation for a pre-defined math file.
 
     Only the changes made relative to the base math will be shown.
 
     Args:
-        base_model (calliope.Model): Calliope model with only the base math applied.
+        base_documentation (MathDocumentation): model documentation with only the base math applied.
         override (str): Name of override to load from the list available in the model config.
+
+    Returns:
+        MathDocumentation: model math documentation with latex backend.
     """
     model = calliope.Model(model_definition=MODEL_PATH, scenario=override)
+    model.build()
 
     full_del = []
     expr_del = []
-    for component_group, component_group_dict in model.math.items():
+    for component_group, component_group_dict in model.applied_math.data.items():
         for name, component_dict in component_group_dict.items():
-            if name in base_model.math[component_group]:
+            if name in base_documentation.math.data[component_group]:
                 if not component_dict.get("active", True):
                     expr_del.append(name)
                     component_dict["description"] = "|REMOVED|"
                     component_dict["active"] = True
-                elif base_model.math[component_group].get(name, {}) != component_dict:
+                elif (
+                    base_documentation.math.data[component_group].get(name, {})
+                    != component_dict
+                ):
                     _add_to_description(component_dict, "|UPDATED|")
                 else:
                     full_del.append(name)
             else:
                 _add_to_description(component_dict, "|NEW|")
 
-    model.math_documentation.build()
+    math_documentation = MathDocumentation(model)
     for key in expr_del:
-        model.math_documentation._instance._dataset[key].attrs["math_string"] = ""
+        math_documentation.backend._dataset[key].attrs["math_string"] = ""
     for key in full_del:
-        del model.math_documentation._instance._dataset[key]
-    for var in model.math_documentation._instance._dataset.values():
+        del math_documentation.backend._dataset[key]
+    for var in math_documentation.backend._dataset.values():
         var.attrs["references"] = var.attrs["references"].intersection(
-            model.math_documentation._instance._dataset.keys()
+            math_documentation.backend._dataset.keys()
         )
         var.attrs["references"] = var.attrs["references"].difference(expr_del)
 
-    logger.info(
-        model.math_documentation._instance._dataset["carrier_in"].attrs["references"]
-    )
+    logger.info(math_documentation.backend._dataset["carrier_in"].attrs["references"])
 
-    return model
+    return math_documentation
 
 
 def _add_to_description(component_dict: dict, update_string: str) -> None:
diff --git a/docs/pre_defined_math/index.md b/docs/pre_defined_math/index.md
index 6a460298..f919fed7 100644
--- a/docs/pre_defined_math/index.md
+++ b/docs/pre_defined_math/index.md
@@ -3,7 +3,7 @@
 As of Calliope version 0.7, the math used to build optimisation problems is stored in YAML files.
 The pre-defined math is a re-implementation of the formerly hardcoded math formulation in this YAML format.
 
-The base math is _always_ applied to your model when you `build` the optimisation problem.
+The pre-defined math for your chosen run [mode](../creating/config.md#configbuildmode) is _always_ applied to your model when you `build` the optimisation problem.
 We have also pre-defined some additional math, which you can _optionally_ load into your model.
 For instance, the [inter-cluster storage][inter-cluster-storage-math] math allows you to track storage levels in technologies more accurately when you are using timeseries clustering in your model.
 
@@ -11,17 +11,12 @@ To load optional, pre-defined math on top of the base math, you can reference it
 
 ```yaml
 config:
-  init:
+  build:
     add_math: [storage_inter_cluster]
 ```
 
-When solving the model in a run mode other than `plan`, some pre-defined additional math will be applied automatically from a file of the same name (e.g., `spores` mode math is stored in [math/spores.yaml](https://github.com/calliope-project/calliope/blob/main/src/calliope/math/spores.yaml)).
-
-!!! note
-
-    Additional math is applied in the order it appears in the `#!yaml config.init.add_math` list.
-    By default, any run mode math will be applied as the final step.
-    If you want to apply your own math *after* the run mode math, you should add the name of the run mode explicitly to the `#!yaml config.init.add_math` list, e.g., `#!yaml config.init.add_math: [operate, user_defined_math.yaml]`.
+If you are running in the `plan` run mode, this will first apply all the [`plan`][base-math] pre-defined math, then the [`storage_inter_cluster`][inter-cluster-storage-math] pre-defined math.
+All pre-defined math YAML files can be found in [`math` directory of the Calliope source code](https://github.com/calliope-project/calliope/blob/main/src/calliope/math/storage_inter_cluster.yaml).
 
 If you want to introduce new constraints, decision variables, or objectives, you can do so as part of the collection of YAML files describing your model.
 See the [user-defined math](../user_defined_math/index.md) section for an in-depth guide to applying your own math.
diff --git a/docs/user_defined_math/components.md b/docs/user_defined_math/components.md
index b180d347..17966728 100644
--- a/docs/user_defined_math/components.md
+++ b/docs/user_defined_math/components.md
@@ -12,7 +12,7 @@ A decision variable in Calliope math looks like this:
 
 ```yaml
 variables:
---8<-- "src/calliope/math/base.yaml:variable"
+--8<-- "src/calliope/math/plan.yaml:variable"
 ```
 
 1. It needs a unique name (`storage_cap` in the example above).
@@ -48,7 +48,7 @@ To not clutter the objective function with all combinations of variables and par
 
 ```yaml
 global_expressions:
---8<-- "src/calliope/math/base.yaml:expression"
+--8<-- "src/calliope/math/plan.yaml:expression"
 ```
 
 Global expressions are by no means necessary to include, but can make more complex linear expressions easier to keep track of and can reduce post-processing requirements.
@@ -74,7 +74,7 @@ Here is an example:
 
 ```yaml
 constraints:
---8<-- "src/calliope/math/base.yaml:constraint"
+--8<-- "src/calliope/math/plan.yaml:constraint"
 ```
 
 1. It needs a unique name (`set_storage_initial` in the above example).
@@ -138,7 +138,7 @@ With your constrained decision variables and a global expression that binds thes
 
 ```yaml
 objectives:
---8<-- "src/calliope/math/base.yaml:objective"
+--8<-- "src/calliope/math/plan.yaml:objective"
 ```
 
 1. It needs a unique name.
diff --git a/docs/user_defined_math/customise.md b/docs/user_defined_math/customise.md
index bdc0c048..19131ba8 100644
--- a/docs/user_defined_math/customise.md
+++ b/docs/user_defined_math/customise.md
@@ -4,8 +4,8 @@ Once you understand the [math components](components.md) and the [formulation sy
 
 You can find examples of additional math that we have put together in our [math example gallery](examples/index.md).
 
-Whenever you introduce your own math, it will be applied on top of the [base math][base-math].
-Therefore, you can include base math overrides as well as add new math.
+Whenever you introduce your own math, it will be applied on top of the pre-defined math for your chosen run [mode](../creating/config.md#configbuildmode).
+Therefore, you can override the pre-defined math as well as add new math.
 For example, you may want to introduce a timeseries parameter to the pre-defined `storage_max` constraint to limit maximum storage capacity on a per-timestep basis:
 
 ```yaml
@@ -16,11 +16,16 @@ storage_max:
 
 The other elements of the `storage_max` constraints have not changed (`foreach`, `where`, ...), so we do not need to define them again when adding our own twist on the pre-defined math.
 
-When defining your model, you can reference any number of YAML files containing the math you want to add in `config.init`. The paths are relative to your main model configuration file:
+!!! note
+
+    If you prefer to start from scratch with your math, you can ask Calliope to _not_ load the pre-defined math for your chosen run mode by setting `#!yaml config.build.ignore_mode_math: true`.
+
+When defining your model, you can reference any number of YAML files containing the math you want to add in `config.build`.
+The paths are relative to your main model configuration file:
 
 ```yaml
 config:
-  init:
+  build:
     add_math: [my_new_math_1.yaml, my_new_math_2.yaml]
 ```
 
@@ -28,10 +33,22 @@ You can also define a mixture of your own math and the [pre-defined math](../pre
 
 ```yaml
 config:
-  init:
+  build:
     add_math: [my_new_math_1.yaml, storage_inter_cluster, my_new_math_2.md]
 ```
 
+Finally, when working in an interactive Python session, you can add math as a dictionary at build time:
+
+```python
+model.build(add_math_dict={...})
+```
+
+This will be applied after the pre-defined mode math and any math from file listed in `config.build.add_math`.
+
+!!! note
+
+    When working in an interactive Python session, you can view the final math dictionary that has been applied to build the optimisation problem by inspecting `model.applied_math` after a successful call to `model.build()`.
+
 ## Adding your parameters to the YAML schema
 
 Our YAML schemas are used to validate user inputs.
@@ -90,9 +107,13 @@ You can write your model's mathematical formulation to view it in a rich-text fo
 To write a LaTeX, reStructuredText, or Markdown file that includes only the math valid for your model:
 
 ```python
+from calliope.postprocess.math_documentation import MathDocumentation
+
 model = calliope.Model("path/to/model.yaml")
-model.math_documentation.build(include="valid")
-model.math_documentation.write(filename="path/to/output/file.[tex|rst|md]")
+model.build()
+
+math_documentation = MathDocumentation(model, include="valid")
+math_documentation.write(filename="path/to/output/file.[tex|rst|md]")
 ```
 
 You can then convert this to a PDF or HTML page using your renderer of choice.
@@ -102,5 +123,5 @@ We recommend you only use HTML as the equations can become too long for a PDF pa
 
     You can add interactive elements to your documentation, if you are planning to host them online using MKDocs.
     This includes tabs to flip between rich-text math and the input YAML snippet, and dropdown lists for math component cross-references.
-    Just set the `mkdocs_features` argument to `True` in `model.math_documentation.write`.
+    Just set the `mkdocs_features` argument to `True` in `math_documentation.write`.
     We use this functionality in our [pre-defined math](../pre_defined_math/index.md).
diff --git a/docs/user_defined_math/syntax.md b/docs/user_defined_math/syntax.md
index e8bc6d38..b2266d23 100644
--- a/docs/user_defined_math/syntax.md
+++ b/docs/user_defined_math/syntax.md
@@ -107,7 +107,7 @@ If you are defining a `constraint`, then you also need to define a comparison op
 You do not need to define the sets of math components in expressions, unless you are actively "slicing" them.
 Behind the scenes, we will make sure that every relevant element of the defined `foreach` sets are matched together when applying the expression (we [merge the underlying xarray DataArrays](https://docs.xarray.dev/en/stable/user-guide/combining.html)).
 Slicing math components involves appending the component with square brackets that contain the slices, e.g. `flow_out[carriers=electricity, nodes=[A, B]]` will slice the `flow_out` decision variable to focus on `electricity` in its `carriers` dimension and only has two nodes (`A` and `B`) on its `nodes` dimension.
-To find out what dimensions you can slice a component on, see your input data (`model.inputs`) for parameters and the definition for decision variables in your loaded math dictionary (`model.math.variables`).
+To find out what dimensions you can slice a component on, see your input data (`model.inputs`) for parameters and the definition for decision variables in your math dictionary.
 
 ## Helper functions
 
diff --git a/src/calliope/attrdict.py b/src/calliope/attrdict.py
index 819657a2..bd94df7b 100644
--- a/src/calliope/attrdict.py
+++ b/src/calliope/attrdict.py
@@ -113,6 +113,7 @@ def _resolve_imports(
         loaded: Self,
         resolve_imports: bool | str,
         base_path: str | Path | None = None,
+        allow_override: bool = False,
     ) -> Self:
         if (
             isinstance(resolve_imports, bool)
@@ -137,7 +138,7 @@ def _resolve_imports(
             path = relative_path(base_path, k)
             imported = cls.from_yaml(path)
             # loaded is added to imported (i.e. it takes precedence)
-            imported.union(loaded_dict)
+            imported.union(loaded_dict, allow_override=allow_override)
             loaded_dict = imported
         # 'import' key itself is no longer needed
         loaded_dict.del_key("import")
@@ -151,7 +152,10 @@ def _resolve_imports(
 
     @classmethod
     def from_yaml(
-        cls, filename: str | Path, resolve_imports: bool | str = True
+        cls,
+        filename: str | Path,
+        resolve_imports: bool | str = True,
+        allow_override: bool = False,
     ) -> Self:
         """Returns an AttrDict initialized from the given path or file path.
 
@@ -168,39 +172,54 @@ def from_yaml(
             filename (str | Path): YAML file.
             resolve_imports (bool | str, optional): top-level `import:` solving option.
                 Defaults to True.
+            allow_override (bool, optional): whether or not to allow overrides of already defined keys.
+                Defaults to False.
 
         Returns:
             Self: constructed AttrDict
         """
         filename = Path(filename)
         loaded = cls(_yaml_load(filename.read_text(encoding="utf-8")))
-        loaded = cls._resolve_imports(loaded, resolve_imports, filename)
+        loaded = cls._resolve_imports(
+            loaded, resolve_imports, filename, allow_override=allow_override
+        )
         return loaded
 
     @classmethod
-    def from_yaml_string(cls, string: str, resolve_imports: bool | str = True) -> Self:
+    def from_yaml_string(
+        cls,
+        string: str,
+        resolve_imports: bool | str = True,
+        allow_override: bool = False,
+    ) -> Self:
         """Returns an AttrDict initialized from the given string.
 
         Input string must be valid YAML.
 
+        If `resolve_imports` is True, top-level `import:` statements
+        are resolved recursively.
+        If `resolve_imports` is False, top-level `import:` statements
+        are treated like any other key and not further processed.
+        If `resolve_imports` is a string, such as `foobar`, import
+        statements underneath that key are resolved, i.e. `foobar.import:`.
+        When resolving import statements, anything defined locally
+        overrides definitions in the imported file.
+
         Args:
             string (str): Valid YAML string.
-            resolve_imports (bool | str, optional):
-                If ``resolve_imports`` is True, top-level ``import:`` statements
-                are resolved recursively.
-                If ``resolve_imports is False, top-level ``import:`` statements
-                are treated like any other key and not further processed.
-                If ``resolve_imports`` is a string, such as ``foobar``, import
-                statements underneath that key are resolved, i.e. ``foobar.import:``.
-                When resolving import statements, anything defined locally
-                overrides definitions in the imported file.
+            resolve_imports (bool | str, optional): top-level `import:` solving option.
+                Defaults to True.
+            allow_override (bool, optional): whether or not to allow overrides of already defined keys.
+                Defaults to False.
 
         Returns:
             calliope.AttrDict:
 
         """
         loaded = cls(_yaml_load(string))
-        loaded = cls._resolve_imports(loaded, resolve_imports)
+        loaded = cls._resolve_imports(
+            loaded, resolve_imports, allow_override=allow_override
+        )
         return loaded
 
     def set_key(self, key, value):
diff --git a/src/calliope/backend/__init__.py b/src/calliope/backend/__init__.py
index 6b715743..d37395d8 100644
--- a/src/calliope/backend/__init__.py
+++ b/src/calliope/backend/__init__.py
@@ -5,23 +5,28 @@
 import xarray as xr
 
 from calliope.backend.gurobi_backend_model import GurobiBackendModel
-from calliope.backend.latex_backend_model import MathDocumentation
+from calliope.backend.latex_backend_model import (
+    ALLOWED_MATH_FILE_FORMATS,
+    LatexBackendModel,
+)
 from calliope.backend.parsing import ParsedBackendComponent
 from calliope.backend.pyomo_backend_model import PyomoBackendModel
 from calliope.exceptions import BackendError
-
-MODEL_BACKENDS = ("pyomo",)
+from calliope.preprocess import CalliopeMath
 
 if TYPE_CHECKING:
     from calliope.backend.backend_model import BackendModel
 
 
-def get_model_backend(name: str, data: xr.Dataset, **kwargs) -> "BackendModel":
+def get_model_backend(
+    name: str, data: xr.Dataset, math: CalliopeMath, **kwargs
+) -> "BackendModel":
     """Assign a backend using the given configuration.
 
     Args:
         name (str): name of the backend to use.
         data (Dataset): model data for the backend.
+        math (CalliopeMath): Calliope math.
         **kwargs: backend keyword arguments corresponding to model.config.build.
 
     Raises:
@@ -32,8 +37,8 @@ def get_model_backend(name: str, data: xr.Dataset, **kwargs) -> "BackendModel":
     """
     match name:
         case "pyomo":
-            return PyomoBackendModel(data, **kwargs)
+            return PyomoBackendModel(data, math, **kwargs)
         case "gurobi":
-            return GurobiBackendModel(data, **kwargs)
+            return GurobiBackendModel(data, math, **kwargs)
         case _:
             raise BackendError(f"Incorrect backend '{name}' requested.")
diff --git a/src/calliope/backend/backend_model.py b/src/calliope/backend/backend_model.py
index 26fd70d5..c52d74ab 100644
--- a/src/calliope/backend/backend_model.py
+++ b/src/calliope/backend/backend_model.py
@@ -4,7 +4,6 @@
 
 from __future__ import annotations
 
-import importlib
 import logging
 import time
 import typing
@@ -32,12 +31,11 @@
 from calliope.backend import helper_functions, parsing
 from calliope.exceptions import warn as model_warn
 from calliope.io import load_config
+from calliope.preprocess.model_math import ORDERED_COMPONENTS_T, CalliopeMath
 from calliope.util.schema import (
-    MATH_SCHEMA,
     MODEL_SCHEMA,
     extract_from_schema,
     update_then_validate_config,
-    validate_dict,
 )
 
 if TYPE_CHECKING:
@@ -46,14 +44,8 @@
 from calliope.exceptions import BackendError
 
 T = TypeVar("T")
-_COMPONENTS_T = Literal[
-    "parameters",
-    "variables",
-    "global_expressions",
-    "constraints",
-    "piecewise_constraints",
-    "objectives",
-]
+ALL_COMPONENTS_T = Literal["parameters", ORDERED_COMPONENTS_T]
+
 
 LOGGER = logging.getLogger(__name__)
 
@@ -61,7 +53,7 @@
 class BackendModelGenerator(ABC):
     """Helper class for backends."""
 
-    _VALID_COMPONENTS: tuple[_COMPONENTS_T, ...] = typing.get_args(_COMPONENTS_T)
+    LID_COMPONENTS: tuple[ALL_COMPONENTS_T, ...] = typing.get_args(ALL_COMPONENTS_T)
     _COMPONENT_ATTR_METADATA = [
         "description",
         "unit",
@@ -77,11 +69,12 @@ class BackendModelGenerator(ABC):
     _PARAM_UNITS = extract_from_schema(MODEL_SCHEMA, "x-unit")
     _PARAM_TYPE = extract_from_schema(MODEL_SCHEMA, "x-type")
 
-    def __init__(self, inputs: xr.Dataset, **kwargs):
+    def __init__(self, inputs: xr.Dataset, math: CalliopeMath, **kwargs):
         """Abstract base class to build a representation of the optimisation problem.
 
         Args:
             inputs (xr.Dataset): Calliope model data.
+            math (CalliopeMath): Calliope math.
             **kwargs (Any): build configuration overrides.
         """
         self._dataset = xr.Dataset()
@@ -90,10 +83,12 @@ def __init__(self, inputs: xr.Dataset, **kwargs):
         self.inputs.attrs["config"]["build"] = update_then_validate_config(
             "build", self.inputs.attrs["config"], **kwargs
         )
-        self._check_inputs()
-
+        self.math: CalliopeMath = deepcopy(math)
         self._solve_logger = logging.getLogger(__name__ + ".<solve>")
 
+        self._check_inputs()
+        self.math.validate()
+
     @abstractmethod
     def add_parameter(
         self, parameter_name: str, parameter_values: xr.DataArray, default: Any = np.nan
@@ -180,7 +175,7 @@ def add_objective(
 
     def log(
         self,
-        component_type: _COMPONENTS_T,
+        component_type: ALL_COMPONENTS_T,
         component_name: str,
         message: str,
         level: Literal["info", "warning", "debug", "error", "critical"] = "debug",
@@ -188,7 +183,7 @@ def log(
         """Log to module-level logger with some prettification of the message.
 
         Args:
-            component_type (_COMPONENTS_T): type of component.
+            component_type (ALL_COMPONENTS_T): type of component.
             component_name (str): name of the component.
             message (str): message to log.
             level (Literal["info", "warning", "debug", "error", "critical"], optional): log level. Defaults to "debug".
@@ -223,20 +218,39 @@ def _check_inputs(self):
             check_results["warn"], check_results["fail"]
         )
 
-    def add_all_math(self):
-        """Parse and all the math stored in the input data."""
-        self._add_run_mode_math()
+    def _validate_math_string_parsing(self) -> None:
+        """Validate that `expression` and `where` strings of the math dictionary can be successfully parsed.
+
+        NOTE: strings are not checked for evaluation validity.
+        Evaluation issues will be raised only on adding a component to the backend.
+        """
+        validation_errors: dict = dict()
+        for component_group in typing.get_args(ORDERED_COMPONENTS_T):
+            for name, dict_ in self.math.data[component_group].items():
+                parsed = parsing.ParsedBackendComponent(component_group, name, dict_)
+                parsed.parse_top_level_where(errors="ignore")
+                parsed.parse_equations(self.valid_component_names, errors="ignore")
+                if not parsed._is_valid:
+                    validation_errors[f"{component_group}:{name}"] = parsed._errors
+
+        if validation_errors:
+            exceptions.print_warnings_and_raise_errors(
+                during="math string parsing (marker indicates where parsing stopped, but may not point to the root cause of the issue)",
+                errors=validation_errors,
+            )
+
+        LOGGER.info("Optimisation Model | Validated math strings.")
+
+    def add_optimisation_components(self) -> None:
+        """Parse math and inputs and set optimisation problem."""
         # The order of adding components matters!
         # 1. Variables, 2. Global Expressions, 3. Constraints, 4. Objectives
-        for components in [
-            "variables",
-            "global_expressions",
-            "constraints",
-            "piecewise_constraints",
-            "objectives",
-        ]:
+        self._add_all_inputs_as_parameters()
+        if self.inputs.attrs["config"]["build"]["pre_validate_math_strings"]:
+            self._validate_math_string_parsing()
+        for components in typing.get_args(ORDERED_COMPONENTS_T):
             component = components.removesuffix("s")
-            for name, dict_ in self.inputs.math[components].items():
+            for name, dict_ in self.math.data[components].items():
                 start = time.time()
                 getattr(self, f"add_{component}")(name, dict_)
                 end = time.time() - start
@@ -245,38 +259,12 @@ def add_all_math(self):
                 )
             LOGGER.info(f"Optimisation Model | {components} | Generated.")
 
-    def _add_run_mode_math(self) -> None:
-        """If not given in the add_math list, override model math with run mode math."""
-        # FIXME: available modes should not be hardcoded here. They should come from a YAML schema.
-        mode = self.inputs.attrs["config"].build.mode
-        add_math = self.inputs.attrs["applied_additional_math"]
-        not_run_mode = {"plan", "operate", "spores"}.difference([mode])
-        run_mode_mismatch = not_run_mode.intersection(add_math)
-        if run_mode_mismatch:
-            exceptions.warn(
-                f"Running in {mode} mode, but run mode(s) {run_mode_mismatch} "
-                "math being loaded from file via the model configuration"
-            )
-
-        if mode != "plan" and mode not in add_math:
-            LOGGER.debug(f"Updating math formulation with {mode} mode math.")
-            filepath = importlib.resources.files("calliope") / "math" / f"{mode}.yaml"
-            self.inputs.math.union(AttrDict.from_yaml(filepath), allow_override=True)
-
-        validate_dict(self.inputs.math, MATH_SCHEMA, "math")
-
     def _add_component(
         self,
         name: str,
         component_dict: Tp,
         component_setter: Callable,
-        component_type: Literal[
-            "variables",
-            "global_expressions",
-            "constraints",
-            "piecewise_constraints",
-            "objectives",
-        ],
+        component_type: ORDERED_COMPONENTS_T,
         break_early: bool = True,
     ) -> parsing.ParsedBackendComponent | None:
         """Generalised function to add a optimisation problem component array to the model.
@@ -286,7 +274,7 @@ def _add_component(
                 this name must be available in the input math provided on initialising the class.
             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"]):
+            component_type (Literal["variables", "global_expressions", "constraints", "piecewise_constraints", "objectives"]):
                 type of the added component.
             break_early (bool, optional): break if the component is not active. Defaults to True.
 
@@ -299,8 +287,8 @@ def _add_component(
         """
         references: set[str] = set()
 
-        if name not in self.inputs.math.get(component_type, {}):
-            self.inputs.math.set_key(f"{component_type}.name", component_dict)
+        if name not in self.math.data[component_type]:
+            self.math.add(AttrDict({f"{component_type}.{name}": component_dict}))
 
         if break_early and not component_dict.get("active", True):
             self.log(
@@ -370,7 +358,7 @@ def _add_component(
         return parsed_component
 
     @abstractmethod
-    def delete_component(self, key: str, component_type: _COMPONENTS_T) -> None:
+    def delete_component(self, key: str, component_type: ALL_COMPONENTS_T) -> None:
         """Delete a list object from the backend model object.
 
         Args:
@@ -379,7 +367,7 @@ def delete_component(self, key: str, component_type: _COMPONENTS_T) -> None:
         """
 
     @abstractmethod
-    def _create_obj_list(self, key: str, component_type: _COMPONENTS_T) -> None:
+    def _create_obj_list(self, key: str, component_type: ALL_COMPONENTS_T) -> None:
         """Attach an empty list object to the backend model object.
 
         The attachment may be a backend-specific subclass of a standard list object.
@@ -432,7 +420,7 @@ def _add_to_dataset(
         self,
         name: str,
         da: xr.DataArray,
-        obj_type: _COMPONENTS_T,
+        obj_type: ALL_COMPONENTS_T,
         unparsed_dict: parsing.UNPARSED_DICTS | dict,
         references: set | None = None,
     ):
@@ -441,7 +429,7 @@ def _add_to_dataset(
         Args:
             name (str): Name of entry in dataset.
             da (xr.DataArray): Data to add.
-            obj_type (_COMPONENTS_T): Type of backend objects in the array.
+            obj_type (ALL_COMPONENTS_T): Type of backend objects in the array.
             unparsed_dict (parsing.UNPARSED_DICTS | dict):
                 Dictionary describing the object being added, from which descriptor
                 attributes will be extracted and added to the array attributes.
@@ -537,7 +525,7 @@ def _apply_func(
             da = tuple(arr.fillna(np.nan) for arr in da)
         return da
 
-    def _raise_error_on_preexistence(self, key: str, obj_type: _COMPONENTS_T):
+    def _raise_error_on_preexistence(self, key: str, obj_type: ALL_COMPONENTS_T):
         """Detect if preexistance errors are present the dataset.
 
         We do not allow any overlap of backend object names since they all have to
@@ -546,7 +534,7 @@ def _raise_error_on_preexistence(self, key: str, obj_type: _COMPONENTS_T):
 
         Args:
             key (str): Backend object name
-            obj_type (Literal["variables", "constraints", "objectives", "parameters", "expressions"]): Object type.
+            obj_type (ALL_COMPONENTS_T): Object type.
 
         Raises:
             BackendError: if `key` already exists in the backend model
@@ -609,7 +597,7 @@ def _filter(val):
         in_math = set(
             name
             for component in ["variables", "global_expressions"]
-            for name in self.inputs.math[component].keys()
+            for name in self.math.data[component]
         )
         return in_data.union(in_math)
 
@@ -617,15 +605,18 @@ def _filter(val):
 class BackendModel(BackendModelGenerator, Generic[T]):
     """Calliope's backend model functionality."""
 
-    def __init__(self, inputs: xr.Dataset, instance: T, **kwargs) -> None:
+    def __init__(
+        self, inputs: xr.Dataset, math: CalliopeMath, instance: T, **kwargs
+    ) -> None:
         """Abstract base class to build backend models that interface with solvers.
 
         Args:
             inputs (xr.Dataset): Calliope model data.
+            math (CalliopeMath): Calliope math.
             instance (T): Interface model instance.
             **kwargs: build configuration overrides.
         """
-        super().__init__(inputs, **kwargs)
+        super().__init__(inputs, math, **kwargs)
         self._instance = instance
         self.shadow_prices: ShadowPrices
         self._has_verbose_strings: bool = False
@@ -1028,19 +1019,13 @@ def _rebuild_references(self, references: set[str]) -> None:
         Args:
             references (set[str]): names of optimisation problem components.
         """
-        ordered_components = [
-            "variables",
-            "global_expressions",
-            "constraints",
-            "objectives",
-        ]
-        for component in ordered_components:
+        for component in typing.get_args(ORDERED_COMPONENTS_T):
             # Rebuild references in the order they are found in the backend dataset
             # which should correspond to the order they were added to the optimisation problem.
             refs = [k for k in getattr(self, component).data_vars if k in references]
             for ref in refs:
                 self.delete_component(ref, component)
-                dict_ = self.inputs.attrs["math"][component][ref]
+                dict_ = self.math.data[component][ref]
                 getattr(self, "add_" + component.removesuffix("s"))(ref, dict_)
 
     def _get_component(self, name: str, component_group: str) -> xr.DataArray:
diff --git a/src/calliope/backend/expression_parser.py b/src/calliope/backend/expression_parser.py
index 49139942..0bd4c676 100644
--- a/src/calliope/backend/expression_parser.py
+++ b/src/calliope/backend/expression_parser.py
@@ -37,6 +37,7 @@
 from typing import TYPE_CHECKING, Any, Literal, overload
 
 import numpy as np
+import pandas as pd
 import pyparsing as pp
 import xarray as xr
 from typing_extensions import NotRequired, TypedDict, Unpack
@@ -788,7 +789,7 @@ def as_array(self) -> xr.DataArray:  # noqa: D102, override
                 evaluated = backend_interface._dataset[self.name]
             except KeyError:
                 evaluated = xr.DataArray(self.name, attrs={"obj_type": "string"})
-        if "default" in evaluated.attrs:
+        if "default" in evaluated.attrs and pd.notna(evaluated.attrs["default"]):
             evaluated = evaluated.fillna(evaluated.attrs["default"])
 
         self.eval_attrs["references"].add(self.name)
diff --git a/src/calliope/backend/gurobi_backend_model.py b/src/calliope/backend/gurobi_backend_model.py
index e5f8096d..2d2e0a48 100644
--- a/src/calliope/backend/gurobi_backend_model.py
+++ b/src/calliope/backend/gurobi_backend_model.py
@@ -17,6 +17,7 @@
 from calliope.backend import backend_model, parsing
 from calliope.exceptions import BackendError, BackendWarning
 from calliope.exceptions import warn as model_warn
+from calliope.preprocess import CalliopeMath
 
 if importlib.util.find_spec("gurobipy") is not None:
     import gurobipy
@@ -40,23 +41,22 @@
 class GurobiBackendModel(backend_model.BackendModel):
     """gurobipy-specific backend functionality."""
 
-    def __init__(self, inputs: xr.Dataset, **kwargs) -> None:
+    def __init__(self, inputs: xr.Dataset, math: CalliopeMath, **kwargs) -> None:
         """Gurobi solver interface class.
 
         Args:
             inputs (xr.Dataset): Calliope model data.
+            math (CalliopeMath): Calliope math.
             **kwargs: passed directly to the solver.
         """
         if importlib.util.find_spec("gurobipy") is None:
             raise ImportError(
                 "Install the `gurobipy` package to build the optimisation problem with the Gurobi backend."
             )
-        super().__init__(inputs, gurobipy.Model(), **kwargs)
+        super().__init__(inputs, math, gurobipy.Model(), **kwargs)
         self._instance: gurobipy.Model
         self.shadow_prices = GurobiShadowPrices(self)
 
-        self._add_all_inputs_as_parameters()
-
     def add_parameter(  # noqa: D102, override
         self, parameter_name: str, parameter_values: xr.DataArray, default: Any = np.nan
     ) -> None:
@@ -275,7 +275,7 @@ def _solve(
 
     def verbose_strings(self) -> None:  # noqa: D102, override
         def __renamer(val, *idx, name: str, attr: str):
-            if pd.notnull(val):
+            if pd.notna(val):
                 new_obj_name = f"{name}[{', '.join(idx)}]"
                 setattr(val, attr, new_obj_name)
 
@@ -389,7 +389,7 @@ def update_variable_bounds(  # noqa: D102, override
                 )
                 continue
 
-            existing_bound_param = self.inputs.attrs["math"].get_key(
+            existing_bound_param = self.math.data.get_key(
                 f"variables.{name}.bounds.{bound_name}", None
             )
             if existing_bound_param in self.parameters:
diff --git a/src/calliope/backend/latex_backend_model.py b/src/calliope/backend/latex_backend_model.py
index c1c86243..c33229b0 100644
--- a/src/calliope/backend/latex_backend_model.py
+++ b/src/calliope/backend/latex_backend_model.py
@@ -5,9 +5,7 @@
 import logging
 import re
 import textwrap
-import typing
-from pathlib import Path
-from typing import Any, Literal, overload
+from typing import Any, Literal
 
 import jinja2
 import numpy as np
@@ -16,120 +14,12 @@
 
 from calliope.backend import backend_model, parsing
 from calliope.exceptions import ModelError
+from calliope.preprocess import CalliopeMath
 
-_ALLOWED_MATH_FILE_FORMATS = Literal["tex", "rst", "md"]
-
-
+ALLOWED_MATH_FILE_FORMATS = Literal["tex", "rst", "md"]
 LOGGER = logging.getLogger(__name__)
 
 
-class MathDocumentation:
-    """For math documentation."""
-
-    def __init__(self) -> None:
-        """Math documentation builder/writer.
-
-        Args:
-            backend_builder (Callable):
-                Method to generate all optimisation problem components on a calliope.backend_model.BackendModel object.
-        """
-        self._inputs: xr.Dataset
-
-    def build(self, include: Literal["all", "valid"] = "all", **kwargs) -> None:
-        """Build string representations of the mathematical formulation using LaTeX math notation, ready to be written with `write`.
-
-        Args:
-            include (Literal["all", "valid"], optional):
-                Defines whether to include all possible math equations ("all") or only
-                those for which at least one index item in the "where" string is valid
-                ("valid"). Defaults to "all".
-            **kwargs: kwargs for the LaTeX backend.
-        """
-        backend = LatexBackendModel(self._inputs, include=include, **kwargs)
-        backend.add_all_math()
-
-        self._instance = backend
-
-    @property
-    def inputs(self):
-        """Getter for backend inputs."""
-        return self._inputs
-
-    @inputs.setter
-    def inputs(self, val: xr.Dataset):
-        """Setter for backend inputs."""
-        self._inputs = val
-
-    # Expecting string if not giving filename.
-    @overload
-    def write(
-        self,
-        filename: Literal[None] = None,
-        mkdocs_features: bool = False,
-        format: _ALLOWED_MATH_FILE_FORMATS | None = None,
-    ) -> str: ...
-
-    # Expecting None (and format arg is not needed) if giving filename.
-    @overload
-    def write(self, filename: str | Path, mkdocs_features: bool = False) -> None: ...
-
-    def write(
-        self,
-        filename: str | Path | None = None,
-        mkdocs_features: bool = False,
-        format: _ALLOWED_MATH_FILE_FORMATS | None = None,
-    ) -> str | None:
-        """Write model documentation.
-
-        `build` must be run beforehand.
-
-        Args:
-            filename (str | Path | None, optional):
-                If given, will write the built mathematical formulation to a file with
-                the given extension as the file format. Defaults to None.
-            mkdocs_features (bool, optional):
-                If True and Markdown docs are being generated, then:
-                - the equations will be on a tab and the original YAML math definition will be on another tab;
-                - the equation cross-references will be given in a drop-down list.
-                Defaults to False.
-            format (_ALLOWED_MATH_FILE_FORMATS | None, optional):
-                Not required if filename is given (as the format will be automatically inferred).
-                Required if expecting a string return from calling this function. The LaTeX math will be embedded in a document of the given format (tex=LaTeX, rst=reStructuredText, md=Markdown).
-                Defaults to None.
-
-        Raises:
-            exceptions.ModelError: Math strings need to be built first (`build`)
-            ValueError: The file format (inferred automatically from `filename` or given by `format`) must be one of ["tex", "rst", "md"].
-
-        Returns:
-            str | None:
-                If `filename` is None, the built mathematical formulation documentation will be returned as a string.
-        """
-        if not hasattr(self, "_instance"):
-            raise ModelError(
-                "Build the documentation (`build`) before trying to write it"
-            )
-
-        if format is None and filename is not None:
-            format = Path(filename).suffix.removeprefix(".")  # type: ignore
-            LOGGER.info(
-                f"Inferring math documentation format from filename as `{format}`."
-            )
-
-        allowed_formats = typing.get_args(_ALLOWED_MATH_FILE_FORMATS)
-        if format is None or format not in allowed_formats:
-            raise ValueError(
-                f"Math documentation format must be one of {allowed_formats}, received `{format}`"
-            )
-        populated_doc = self._instance.generate_math_doc(format, mkdocs_features)
-
-        if filename is None:
-            return populated_doc
-        else:
-            Path(filename).write_text(populated_doc)
-            return None
-
-
 class LatexBackendModel(backend_model.BackendModelGenerator):
     """Calliope's LaTeX backend."""
 
@@ -412,21 +302,24 @@ class LatexBackendModel(backend_model.BackendModelGenerator):
     FORMAT_STRINGS = {"rst": RST_DOC, "tex": TEX_DOC, "md": MD_DOC}
 
     def __init__(
-        self, inputs: xr.Dataset, include: Literal["all", "valid"] = "all", **kwargs
+        self,
+        inputs: xr.Dataset,
+        math: CalliopeMath,
+        include: Literal["all", "valid"] = "all",
+        **kwargs,
     ) -> None:
         """Interface to build a string representation of the mathematical formulation using LaTeX math notation.
 
         Args:
             inputs (xr.Dataset): model data.
+            math (CalliopeMath): Calliope math.
             include (Literal["all", "valid"], optional):
                 Defines whether to include all possible math equations ("all") or only those for which at least one index item in the "where" string is valid ("valid"). Defaults to "all".
             **kwargs: for the backend model generator.
         """
-        super().__init__(inputs, **kwargs)
+        super().__init__(inputs, math, **kwargs)
         self.include = include
 
-        self._add_all_inputs_as_parameters()
-
     def add_parameter(  # noqa: D102, override
         self, parameter_name: str, parameter_values: xr.DataArray, default: Any = np.nan
     ) -> None:
@@ -539,8 +432,6 @@ def add_variable(  # noqa: D102, override
         def _variable_setter(where: xr.DataArray, references: set) -> xr.DataArray:
             return where.where(where)
 
-        domain = domain_dict[variable_dict.get("domain", "real")]
-
         parsed_component = self._add_component(
             name, variable_dict, _variable_setter, "variables", break_early=False
         )
@@ -568,7 +459,6 @@ def add_objective(  # noqa: D102, override
             "minimise": r"\min{}",
             "maximise": r"\max{}",
         }
-
         equation_strings: list = []
 
         def _objective_setter(
@@ -599,7 +489,7 @@ def delete_component(  # noqa: D102, override
             del self._dataset[key]
 
     def generate_math_doc(
-        self, format: _ALLOWED_MATH_FILE_FORMATS = "tex", mkdocs_features: bool = False
+        self, format: ALLOWED_MATH_FILE_FORMATS = "tex", mkdocs_features: bool = False
     ) -> str:
         """Generate the math documentation by embedding LaTeX math in a template.
 
@@ -657,7 +547,7 @@ def generate_math_doc(
             ]
             if getattr(self, objtype).data_vars
         }
-        if not components["parameters"]:
+        if "parameters" in components and not components["parameters"]:
             del components["parameters"]
         return self._render(
             doc_template, mkdocs_features=mkdocs_features, components=components
diff --git a/src/calliope/backend/parsing.py b/src/calliope/backend/parsing.py
index 278f037a..33c9ea47 100644
--- a/src/calliope/backend/parsing.py
+++ b/src/calliope/backend/parsing.py
@@ -850,13 +850,13 @@ def generate_top_level_where_array(
         return where
 
     def raise_caught_errors(self):
-        """If there are any parsing errors, pipe them to the ModelError bullet point list generator."""
+        """Pipe parsing errors to the ModelError bullet point list generator."""
         if not self._is_valid:
             exceptions.print_warnings_and_raise_errors(
                 errors={f"{self.name}": self._errors},
                 during=(
                     "math string parsing (marker indicates where parsing stopped, "
-                    "which might not be the root cause of the issue; sorry...)"
+                    "but may not point to the root cause of the issue)"
                 ),
                 bullet=self._ERR_BULLET,
             )
diff --git a/src/calliope/backend/pyomo_backend_model.py b/src/calliope/backend/pyomo_backend_model.py
index 228b085a..5ba41ba0 100644
--- a/src/calliope/backend/pyomo_backend_model.py
+++ b/src/calliope/backend/pyomo_backend_model.py
@@ -26,11 +26,13 @@
 from pyomo.opt import SolverFactory  # type: ignore
 from pyomo.util.model_size import build_model_size_report  # type: ignore
 
-from calliope.backend import backend_model, parsing
 from calliope.exceptions import BackendError, BackendWarning
 from calliope.exceptions import warn as model_warn
+from calliope.preprocess import CalliopeMath
 from calliope.util.logging import LogWriter
 
+from . import backend_model, parsing
+
 T = TypeVar("T")
 _COMPONENTS_T = Literal[
     "variables",
@@ -56,14 +58,15 @@
 class PyomoBackendModel(backend_model.BackendModel):
     """Pyomo-specific backend functionality."""
 
-    def __init__(self, inputs: xr.Dataset, **kwargs) -> None:
+    def __init__(self, inputs: xr.Dataset, math: CalliopeMath, **kwargs) -> None:
         """Pyomo solver interface class.
 
         Args:
             inputs (xr.Dataset): Calliope model data.
+            math (CalliopeMath): Calliope math.
             **kwargs: passed directly to the solver.
         """
-        super().__init__(inputs, pmo.block(), **kwargs)
+        super().__init__(inputs, math, pmo.block(), **kwargs)
 
         self._instance.parameters = pmo.parameter_dict()
         self._instance.variables = pmo.variable_dict()
@@ -75,8 +78,6 @@ def __init__(self, inputs: xr.Dataset, **kwargs) -> None:
         self._instance.dual = pmo.suffix(direction=pmo.suffix.IMPORT)
         self.shadow_prices = PyomoShadowPrices(self._instance.dual, self)
 
-        self._add_all_inputs_as_parameters()
-
     def add_parameter(  # noqa: D102, override
         self, parameter_name: str, parameter_values: xr.DataArray, default: Any = np.nan
     ) -> None:
@@ -330,7 +331,7 @@ def _solve(  # noqa: D102, override
 
     def verbose_strings(self) -> None:  # noqa: D102, override
         def __renamer(val, *idx):
-            if pd.notnull(val):
+            if pd.notna(val):
                 val.calliope_coords = idx
 
         with self._datetime_as_string(self._dataset):
@@ -460,7 +461,7 @@ def update_variable_bounds(  # noqa: D102, override
                 )
                 continue
 
-            existing_bound_param = self.inputs.attrs["math"].get_key(
+            existing_bound_param = self.math.data.get_key(
                 f"variables.{name}.bounds.{bound_name}", None
             )
             if existing_bound_param in self.parameters:
diff --git a/src/calliope/config/config_schema.yaml b/src/calliope/config/config_schema.yaml
index 463dbe4c..b9ebe627 100644
--- a/src/calliope/config/config_schema.yaml
+++ b/src/calliope/config/config_schema.yaml
@@ -51,17 +51,6 @@ properties:
             type: string
             default: "ISO8601"
             description: Timestamp format of all time series data when read from file. "ISO8601" means "%Y-%m-%d %H:%M:%S".
-          add_math:
-            type: array
-            default: []
-            description: List of references to files which contain additional mathematical formulations to be applied on top of the base math.
-            uniqueItems: true
-            items:
-              type: string
-              description: >
-                If referring to an pre-defined Calliope math file (see documentation for available files), do not append the reference with ".yaml".
-                If referring to your own math file, ensure the file type is given as a suffix (".yaml" or ".yml").
-                Relative paths will be assumed to be relative to the model definition file given when creating a calliope Model (`calliope.Model(model_definition=...)`)
           distance_unit:
             type: string
             default: km
@@ -77,6 +66,23 @@ properties:
           Additional configuration items will be passed onto math string parsing and can therefore be accessed in the `where` strings by `config.[item-name]`, where "[item-name]" is the name of your own configuration item.
         additionalProperties: true
         properties:
+          add_math:
+            type: array
+            default: []
+            description: List of references to files which contain additional mathematical formulations to be applied on top of or instead of the base mode math.
+            uniqueItems: true
+            items:
+              type: string
+              description: >
+                If referring to an pre-defined Calliope math file (see documentation for available files), do not append the reference with ".yaml".
+                If referring to your own math file, ensure the file type is given as a suffix (".yaml" or ".yml").
+                Relative paths will be assumed to be relative to the model definition file given when creating a calliope Model (`calliope.Model(model_definition=...)`).
+          ignore_mode_math:
+            type: boolean
+            default: false
+            description: >-
+              If True, do not initialise the mathematical formulation with the pre-defined math for the given run `mode`.
+              This option can be used to completely re-define the Calliope mathematical formulation.
           backend:
             type: string
             default: pyomo
@@ -111,6 +117,12 @@ properties:
             type: boolean
             default: false
             description: If the model already contains `plan` mode results, use those optimal capacities as input parameters to the `operate` mode run.
+          pre_validate_math_strings:
+            type: boolean
+            default: true
+            description: >-
+              If true, the Calliope math definition will be scanned for parsing errors _before_ undertaking the much more expensive operation of building the optimisation problem.
+              You can switch this off (e.g., if you know there are no parsing errors) to reduce overall build time.
 
       solve:
         type: object
diff --git a/src/calliope/example_models/national_scale/data_tables/cluster_days.csv b/src/calliope/example_models/national_scale/data_tables/cluster_days.csv
index 4cd9e78d..587e3dd8 100644
--- a/src/calliope/example_models/national_scale/data_tables/cluster_days.csv
+++ b/src/calliope/example_models/national_scale/data_tables/cluster_days.csv
@@ -1,366 +1,366 @@
-datesteps,cluster
-2005-01-01,2005-12-09
-2005-01-02,2005-12-08
-2005-01-03,2005-12-08
-2005-01-04,2005-12-09
-2005-01-05,2005-12-09
-2005-01-06,2005-12-06
-2005-01-07,2005-12-06
-2005-01-08,2005-12-06
-2005-01-09,2005-12-06
-2005-01-10,2005-12-06
-2005-01-11,2005-12-06
-2005-01-12,2005-12-06
-2005-01-13,2005-12-06
-2005-01-14,2005-12-09
-2005-01-15,2005-12-09
-2005-01-16,2005-12-08
-2005-01-17,2005-12-08
-2005-01-18,2005-12-06
-2005-01-19,2005-12-06
-2005-01-20,2005-12-06
-2005-01-21,2005-12-08
-2005-01-22,2005-12-08
-2005-01-23,2005-12-08
-2005-01-24,2005-12-06
-2005-01-25,2005-12-08
-2005-01-26,2005-12-06
-2005-01-27,2005-12-06
-2005-01-28,2005-12-06
-2005-01-29,2005-12-06
-2005-01-30,2005-12-08
-2005-01-31,2005-02-17
-2005-02-01,2005-12-08
-2005-02-02,2005-12-06
-2005-02-03,2005-02-17
-2005-02-04,2005-12-06
-2005-02-05,2005-12-09
-2005-02-06,2005-12-09
-2005-02-07,2005-02-17
-2005-02-08,2005-12-06
-2005-02-09,2005-12-06
-2005-02-10,2005-12-06
-2005-02-11,2005-12-06
-2005-02-12,2005-12-06
-2005-02-13,2005-12-06
-2005-02-14,2005-12-08
-2005-02-15,2005-02-17
-2005-02-16,2005-12-08
-2005-02-17,2005-12-06
-2005-02-18,2005-12-08
-2005-02-19,2005-12-06
-2005-02-20,2005-12-06
-2005-02-21,2005-02-17
-2005-02-22,2005-12-06
-2005-02-23,2005-12-09
-2005-02-24,2005-02-17
-2005-02-25,2005-09-13
-2005-02-26,2005-04-13
-2005-02-27,2005-12-08
-2005-02-28,2005-12-08
-2005-03-01,2005-12-08
-2005-03-02,2005-12-09
-2005-03-03,2005-02-17
-2005-03-04,2005-09-13
-2005-03-05,2005-12-08
-2005-03-06,2005-02-17
-2005-03-07,2005-04-13
-2005-03-08,2005-04-13
-2005-03-09,2005-04-13
-2005-03-10,2005-04-13
-2005-03-11,2005-04-13
-2005-03-12,2005-04-27
-2005-03-13,2005-04-27
-2005-03-14,2005-12-06
-2005-03-15,2005-09-13
-2005-03-16,2005-09-13
-2005-03-17,2005-09-13
-2005-03-18,2005-09-13
-2005-03-19,2005-04-27
-2005-03-20,2005-04-13
-2005-03-21,2005-04-27
-2005-03-22,2005-04-13
-2005-03-23,2005-04-27
-2005-03-24,2005-09-13
-2005-03-25,2005-09-13
-2005-03-26,2005-04-13
-2005-03-27,2005-04-13
-2005-03-28,2005-12-06
-2005-03-29,2005-04-13
-2005-03-30,2005-04-13
-2005-03-31,2005-04-13
-2005-04-01,2005-04-13
-2005-04-02,2005-12-09
-2005-04-03,2005-12-08
-2005-04-04,2005-09-13
-2005-04-05,2005-09-13
-2005-04-06,2005-09-13
-2005-04-07,2005-04-13
-2005-04-08,2005-12-09
-2005-04-09,2005-08-09
-2005-04-10,2005-08-09
-2005-04-11,2005-08-09
-2005-04-12,2005-08-09
-2005-04-13,2005-08-09
-2005-04-14,2005-04-13
-2005-04-15,2005-04-13
-2005-04-16,2005-09-13
-2005-04-17,2005-09-09
-2005-04-18,2005-05-25
-2005-04-19,2005-08-09
-2005-04-20,2005-08-09
-2005-04-21,2005-09-09
-2005-04-22,2005-09-09
-2005-04-23,2005-09-09
-2005-04-24,2005-09-13
-2005-04-25,2005-08-09
-2005-04-26,2005-08-09
-2005-04-27,2005-08-09
-2005-04-28,2005-08-09
-2005-04-29,2005-09-09
-2005-04-30,2005-09-09
-2005-05-01,2005-09-09
-2005-05-02,2005-09-09
-2005-05-03,2005-09-09
-2005-05-04,2005-08-09
-2005-05-05,2005-08-09
-2005-05-06,2005-09-09
-2005-05-07,2005-08-09
-2005-05-08,2005-09-09
-2005-05-09,2005-09-13
-2005-05-10,2005-08-09
-2005-05-11,2005-04-27
-2005-05-12,2005-04-27
-2005-05-13,2005-04-27
-2005-05-14,2005-09-09
-2005-05-15,2005-08-09
-2005-05-16,2005-04-27
-2005-05-17,2005-08-09
-2005-05-18,2005-08-09
-2005-05-19,2005-08-09
-2005-05-20,2005-09-09
-2005-05-21,2005-05-25
-2005-05-22,2005-08-09
-2005-05-23,2005-09-09
-2005-05-24,2005-08-09
-2005-05-25,2005-09-09
-2005-05-26,2005-09-09
-2005-05-27,2005-09-09
-2005-05-28,2005-09-09
-2005-05-29,2005-04-27
-2005-05-30,2005-04-27
-2005-05-31,2005-12-09
-2005-06-01,2005-09-09
-2005-06-02,2005-09-09
-2005-06-03,2005-08-09
-2005-06-04,2005-08-09
-2005-06-05,2005-08-09
-2005-06-06,2005-08-09
-2005-06-07,2005-08-09
-2005-06-08,2005-08-09
-2005-06-09,2005-08-09
-2005-06-10,2005-08-09
-2005-06-11,2005-09-09
-2005-06-12,2005-08-09
-2005-06-13,2005-09-13
-2005-06-14,2005-08-09
-2005-06-15,2005-05-25
-2005-06-16,2005-08-09
-2005-06-17,2005-08-09
-2005-06-18,2005-08-09
-2005-06-19,2005-08-09
-2005-06-20,2005-08-09
-2005-06-21,2005-08-09
-2005-06-22,2005-09-13
-2005-06-23,2005-04-27
-2005-06-24,2005-08-09
-2005-06-25,2005-08-09
-2005-06-26,2005-08-09
-2005-06-27,2005-04-27
-2005-06-28,2005-08-09
-2005-06-29,2005-09-09
-2005-06-30,2005-08-09
-2005-07-01,2005-08-09
-2005-07-02,2005-08-09
-2005-07-03,2005-08-09
-2005-07-04,2005-09-09
-2005-07-05,2005-09-09
-2005-07-06,2005-08-09
-2005-07-07,2005-08-09
-2005-07-08,2005-08-09
-2005-07-09,2005-08-09
-2005-07-10,2005-08-09
-2005-07-11,2005-08-09
-2005-07-12,2005-08-09
-2005-07-13,2005-08-09
-2005-07-14,2005-08-09
-2005-07-15,2005-09-09
-2005-07-16,2005-08-09
-2005-07-17,2005-08-09
-2005-07-18,2005-08-09
-2005-07-19,2005-08-09
-2005-07-20,2005-08-09
-2005-07-21,2005-08-09
-2005-07-22,2005-08-09
-2005-07-23,2005-08-09
-2005-07-24,2005-08-09
-2005-07-25,2005-08-09
-2005-07-26,2005-08-09
-2005-07-27,2005-09-13
-2005-07-28,2005-09-13
-2005-07-29,2005-08-09
-2005-07-30,2005-08-09
-2005-07-31,2005-08-09
-2005-08-01,2005-08-09
-2005-08-02,2005-09-09
-2005-08-03,2005-08-09
-2005-08-04,2005-08-09
-2005-08-05,2005-09-09
-2005-08-06,2005-09-09
-2005-08-07,2005-09-09
-2005-08-08,2005-09-13
-2005-08-09,2005-09-13
-2005-08-10,2005-04-27
-2005-08-11,2005-05-25
-2005-08-12,2005-08-09
-2005-08-13,2005-08-09
-2005-08-14,2005-08-09
-2005-08-15,2005-08-09
-2005-08-16,2005-09-13
-2005-08-17,2005-08-09
-2005-08-18,2005-08-09
-2005-08-19,2005-04-27
-2005-08-20,2005-08-09
-2005-08-21,2005-08-09
-2005-08-22,2005-08-09
-2005-08-23,2005-08-09
-2005-08-24,2005-08-09
-2005-08-25,2005-09-13
-2005-08-26,2005-08-09
-2005-08-27,2005-08-09
-2005-08-28,2005-05-25
-2005-08-29,2005-09-13
-2005-08-30,2005-08-09
-2005-08-31,2005-09-13
-2005-09-01,2005-08-09
-2005-09-02,2005-08-09
-2005-09-03,2005-05-25
-2005-09-04,2005-09-13
-2005-09-05,2005-09-13
-2005-09-06,2005-09-13
-2005-09-07,2005-04-27
-2005-09-08,2005-09-13
-2005-09-09,2005-05-25
-2005-09-10,2005-08-09
-2005-09-11,2005-09-13
-2005-09-12,2005-09-13
-2005-09-13,2005-09-13
-2005-09-14,2005-08-09
-2005-09-15,2005-09-13
-2005-09-16,2005-04-13
-2005-09-17,2005-04-27
-2005-09-18,2005-09-13
-2005-09-19,2005-09-13
-2005-09-20,2005-09-13
-2005-09-21,2005-09-13
-2005-09-22,2005-09-13
-2005-09-23,2005-12-09
-2005-09-24,2005-09-09
-2005-09-25,2005-09-13
-2005-09-26,2005-05-25
-2005-09-27,2005-05-25
-2005-09-28,2005-05-25
-2005-09-29,2005-09-13
-2005-09-30,2005-05-25
-2005-10-01,2005-05-25
-2005-10-02,2005-04-13
-2005-10-03,2005-09-13
-2005-10-04,2005-09-13
-2005-10-05,2005-09-13
-2005-10-06,2005-09-13
-2005-10-07,2005-09-13
-2005-10-08,2005-04-27
-2005-10-09,2005-04-27
-2005-10-10,2005-12-09
-2005-10-11,2005-12-09
-2005-10-12,2005-04-27
-2005-10-13,2005-04-27
-2005-10-14,2005-12-06
-2005-10-15,2005-09-13
-2005-10-16,2005-12-08
-2005-10-17,2005-02-17
-2005-10-18,2005-04-27
-2005-10-19,2005-09-13
-2005-10-20,2005-04-27
-2005-10-21,2005-05-25
-2005-10-22,2005-04-27
-2005-10-23,2005-09-13
-2005-10-24,2005-09-13
-2005-10-25,2005-09-13
-2005-10-26,2005-09-13
-2005-10-27,2005-09-13
-2005-10-28,2005-04-27
-2005-10-29,2005-09-13
-2005-10-30,2005-12-06
-2005-10-31,2005-09-13
-2005-11-01,2005-04-27
-2005-11-02,2005-12-08
-2005-11-03,2005-12-09
-2005-11-04,2005-12-08
-2005-11-05,2005-12-08
-2005-11-06,2005-12-08
-2005-11-07,2005-12-06
-2005-11-08,2005-09-13
-2005-11-09,2005-09-13
-2005-11-10,2005-04-27
-2005-11-11,2005-04-27
-2005-11-12,2005-04-27
-2005-11-13,2005-12-09
-2005-11-14,2005-04-27
-2005-11-15,2005-04-27
-2005-11-16,2005-12-06
-2005-11-17,2005-12-09
-2005-11-18,2005-12-09
-2005-11-19,2005-12-09
-2005-11-20,2005-12-09
-2005-11-21,2005-12-09
-2005-11-22,2005-12-06
-2005-11-23,2005-12-08
-2005-11-24,2005-12-08
-2005-11-25,2005-02-17
-2005-11-26,2005-12-08
-2005-11-27,2005-12-08
-2005-11-28,2005-12-06
-2005-11-29,2005-12-09
-2005-11-30,2005-12-06
-2005-12-01,2005-12-06
-2005-12-02,2005-12-09
-2005-12-03,2005-12-06
-2005-12-04,2005-12-06
-2005-12-05,2005-12-06
-2005-12-06,2005-12-06
-2005-12-07,2005-12-06
-2005-12-08,2005-12-06
-2005-12-09,2005-12-06
-2005-12-10,2005-12-06
-2005-12-11,2005-12-06
-2005-12-12,2005-12-06
-2005-12-13,2005-12-06
-2005-12-14,2005-12-06
-2005-12-15,2005-12-08
-2005-12-16,2005-12-08
-2005-12-17,2005-12-06
-2005-12-18,2005-12-08
-2005-12-19,2005-12-09
-2005-12-20,2005-12-09
-2005-12-21,2005-12-06
-2005-12-22,2005-12-06
-2005-12-23,2005-12-09
-2005-12-24,2005-12-09
-2005-12-25,2005-12-09
-2005-12-26,2005-12-09
-2005-12-27,2005-12-09
-2005-12-28,2005-12-09
-2005-12-29,2005-12-06
-2005-12-30,2005-12-09
-2005-12-31,2005-12-09
+timesteps,cluster
+2005-01-01,2005-11-04
+2005-01-02,2005-11-04
+2005-01-03,2005-12-07
+2005-01-04,2005-12-07
+2005-01-05,2005-12-07
+2005-01-06,2005-12-07
+2005-01-07,2005-12-07
+2005-01-08,2005-11-01
+2005-01-09,2005-12-07
+2005-01-10,2005-12-07
+2005-01-11,2005-12-07
+2005-01-12,2005-12-07
+2005-01-13,2005-12-07
+2005-01-14,2005-12-07
+2005-01-15,2005-11-04
+2005-01-16,2005-12-07
+2005-01-17,2005-12-07
+2005-01-18,2005-12-07
+2005-01-19,2005-12-07
+2005-01-20,2005-12-07
+2005-01-21,2005-11-04
+2005-01-22,2005-11-04
+2005-01-23,2005-12-07
+2005-01-24,2005-12-07
+2005-01-25,2005-12-07
+2005-01-26,2005-12-07
+2005-01-27,2005-12-07
+2005-01-28,2005-12-07
+2005-01-29,2005-11-01
+2005-01-30,2005-12-07
+2005-01-31,2005-12-07
+2005-02-01,2005-12-07
+2005-02-02,2005-12-07
+2005-02-03,2005-12-07
+2005-02-04,2005-12-07
+2005-02-05,2005-11-04
+2005-02-06,2005-12-07
+2005-02-07,2005-12-07
+2005-02-08,2005-12-07
+2005-02-09,2005-12-07
+2005-02-10,2005-12-07
+2005-02-11,2005-11-01
+2005-02-12,2005-11-01
+2005-02-13,2005-12-07
+2005-02-14,2005-12-07
+2005-02-15,2005-12-07
+2005-02-16,2005-12-07
+2005-02-17,2005-12-07
+2005-02-18,2005-11-04
+2005-02-19,2005-11-01
+2005-02-20,2005-12-07
+2005-02-21,2005-12-07
+2005-02-22,2005-12-07
+2005-02-23,2005-12-07
+2005-02-24,2005-12-07
+2005-02-25,2005-11-01
+2005-02-26,2005-11-04
+2005-02-27,2005-12-07
+2005-02-28,2005-12-07
+2005-03-01,2005-12-07
+2005-03-02,2005-12-07
+2005-03-03,2005-12-07
+2005-03-04,2005-11-01
+2005-03-05,2005-11-04
+2005-03-06,2005-12-07
+2005-03-07,2005-03-09
+2005-03-08,2005-03-09
+2005-03-09,2005-03-09
+2005-03-10,2005-03-09
+2005-03-11,2005-03-09
+2005-03-12,2005-11-01
+2005-03-13,2005-11-01
+2005-03-14,2005-03-09
+2005-03-15,2005-03-09
+2005-03-16,2005-03-09
+2005-03-17,2005-03-09
+2005-03-18,2005-09-19
+2005-03-19,2005-11-01
+2005-03-20,2005-03-09
+2005-03-21,2005-11-01
+2005-03-22,2005-03-09
+2005-03-23,2005-11-01
+2005-03-24,2005-09-19
+2005-03-25,2005-11-01
+2005-03-26,2005-11-04
+2005-03-27,2005-03-09
+2005-03-28,2005-03-09
+2005-03-29,2005-03-09
+2005-03-30,2005-03-09
+2005-03-31,2005-03-09
+2005-04-01,2005-11-04
+2005-04-02,2005-11-04
+2005-04-03,2005-11-04
+2005-04-04,2005-09-19
+2005-04-05,2005-09-19
+2005-04-06,2005-09-19
+2005-04-07,2005-03-09
+2005-04-08,2005-11-04
+2005-04-09,2005-09-19
+2005-04-10,2005-09-19
+2005-04-11,2005-09-19
+2005-04-12,2005-09-19
+2005-04-13,2005-09-19
+2005-04-14,2005-11-04
+2005-04-15,2005-11-04
+2005-04-16,2005-09-19
+2005-04-17,2005-05-02
+2005-04-18,2005-09-28
+2005-04-19,2005-09-19
+2005-04-20,2005-09-19
+2005-04-21,2005-05-02
+2005-04-22,2005-05-02
+2005-04-23,2005-05-02
+2005-04-24,2005-09-19
+2005-04-25,2005-09-19
+2005-04-26,2005-09-19
+2005-04-27,2005-09-19
+2005-04-28,2005-09-19
+2005-04-29,2005-05-02
+2005-04-30,2005-05-02
+2005-05-01,2005-05-02
+2005-05-02,2005-05-02
+2005-05-03,2005-05-02
+2005-05-04,2005-09-19
+2005-05-05,2005-09-19
+2005-05-06,2005-05-02
+2005-05-07,2005-09-19
+2005-05-08,2005-05-02
+2005-05-09,2005-05-13
+2005-05-10,2005-05-13
+2005-05-11,2005-05-13
+2005-05-12,2005-05-13
+2005-05-13,2005-05-13
+2005-05-14,2005-05-02
+2005-05-15,2005-09-19
+2005-05-16,2005-05-13
+2005-05-17,2005-09-19
+2005-05-18,2005-09-19
+2005-05-19,2005-09-19
+2005-05-20,2005-05-02
+2005-05-21,2005-09-28
+2005-05-22,2005-09-19
+2005-05-23,2005-05-02
+2005-05-24,2005-09-19
+2005-05-25,2005-05-02
+2005-05-26,2005-05-02
+2005-05-27,2005-05-02
+2005-05-28,2005-05-02
+2005-05-29,2005-05-13
+2005-05-30,2005-05-13
+2005-05-31,2005-11-04
+2005-06-01,2005-08-17
+2005-06-02,2005-08-17
+2005-06-03,2005-08-17
+2005-06-04,2005-08-17
+2005-06-05,2005-08-17
+2005-06-06,2005-08-17
+2005-06-07,2005-08-17
+2005-06-08,2005-08-17
+2005-06-09,2005-08-17
+2005-06-10,2005-08-17
+2005-06-11,2005-08-17
+2005-06-12,2005-06-13
+2005-06-13,2005-06-13
+2005-06-14,2005-08-17
+2005-06-15,2005-08-17
+2005-06-16,2005-08-17
+2005-06-17,2005-08-17
+2005-06-18,2005-08-17
+2005-06-19,2005-08-17
+2005-06-20,2005-08-17
+2005-06-21,2005-06-13
+2005-06-22,2005-06-13
+2005-06-23,2005-06-13
+2005-06-24,2005-08-17
+2005-06-25,2005-08-17
+2005-06-26,2005-08-17
+2005-06-27,2005-06-13
+2005-06-28,2005-08-17
+2005-06-29,2005-08-17
+2005-06-30,2005-08-17
+2005-07-01,2005-08-17
+2005-07-02,2005-08-17
+2005-07-03,2005-08-17
+2005-07-04,2005-08-17
+2005-07-05,2005-08-17
+2005-07-06,2005-08-17
+2005-07-07,2005-08-17
+2005-07-08,2005-08-17
+2005-07-09,2005-08-17
+2005-07-10,2005-08-17
+2005-07-11,2005-08-17
+2005-07-12,2005-08-17
+2005-07-13,2005-08-17
+2005-07-14,2005-08-17
+2005-07-15,2005-08-17
+2005-07-16,2005-08-17
+2005-07-17,2005-08-17
+2005-07-18,2005-08-17
+2005-07-19,2005-08-17
+2005-07-20,2005-08-17
+2005-07-21,2005-08-17
+2005-07-22,2005-08-17
+2005-07-23,2005-08-17
+2005-07-24,2005-08-17
+2005-07-25,2005-08-17
+2005-07-26,2005-08-17
+2005-07-27,2005-06-13
+2005-07-28,2005-06-13
+2005-07-29,2005-08-17
+2005-07-30,2005-08-17
+2005-07-31,2005-08-17
+2005-08-01,2005-08-17
+2005-08-02,2005-08-17
+2005-08-03,2005-08-17
+2005-08-04,2005-08-17
+2005-08-05,2005-08-17
+2005-08-06,2005-08-17
+2005-08-07,2005-08-17
+2005-08-08,2005-06-13
+2005-08-09,2005-06-13
+2005-08-10,2005-06-13
+2005-08-11,2005-06-13
+2005-08-12,2005-08-17
+2005-08-13,2005-08-17
+2005-08-14,2005-08-17
+2005-08-15,2005-08-17
+2005-08-16,2005-06-13
+2005-08-17,2005-08-17
+2005-08-18,2005-08-17
+2005-08-19,2005-06-13
+2005-08-20,2005-08-17
+2005-08-21,2005-08-17
+2005-08-22,2005-08-17
+2005-08-23,2005-08-17
+2005-08-24,2005-08-17
+2005-08-25,2005-06-13
+2005-08-26,2005-08-17
+2005-08-27,2005-08-17
+2005-08-28,2005-08-17
+2005-08-29,2005-08-17
+2005-08-30,2005-08-17
+2005-08-31,2005-06-13
+2005-09-01,2005-09-19
+2005-09-02,2005-09-19
+2005-09-03,2005-09-28
+2005-09-04,2005-09-19
+2005-09-05,2005-09-19
+2005-09-06,2005-09-19
+2005-09-07,2005-11-01
+2005-09-08,2005-09-19
+2005-09-09,2005-09-28
+2005-09-10,2005-09-19
+2005-09-11,2005-09-19
+2005-09-12,2005-09-19
+2005-09-13,2005-09-19
+2005-09-14,2005-09-19
+2005-09-15,2005-09-19
+2005-09-16,2005-11-04
+2005-09-17,2005-05-13
+2005-09-18,2005-09-19
+2005-09-19,2005-09-19
+2005-09-20,2005-09-19
+2005-09-21,2005-09-19
+2005-09-22,2005-09-19
+2005-09-23,2005-11-04
+2005-09-24,2005-05-02
+2005-09-25,2005-09-19
+2005-09-26,2005-09-28
+2005-09-27,2005-09-28
+2005-09-28,2005-09-28
+2005-09-29,2005-09-19
+2005-09-30,2005-09-28
+2005-10-01,2005-09-28
+2005-10-02,2005-11-04
+2005-10-03,2005-09-19
+2005-10-04,2005-09-19
+2005-10-05,2005-09-19
+2005-10-06,2005-09-19
+2005-10-07,2005-09-19
+2005-10-08,2005-05-13
+2005-10-09,2005-11-01
+2005-10-10,2005-11-04
+2005-10-11,2005-11-04
+2005-10-12,2005-11-01
+2005-10-13,2005-11-01
+2005-10-14,2005-11-01
+2005-10-15,2005-09-19
+2005-10-16,2005-11-04
+2005-10-17,2005-11-04
+2005-10-18,2005-11-01
+2005-10-19,2005-03-09
+2005-10-20,2005-11-01
+2005-10-21,2005-09-28
+2005-10-22,2005-05-13
+2005-10-23,2005-11-01
+2005-10-24,2005-11-01
+2005-10-25,2005-11-01
+2005-10-26,2005-11-01
+2005-10-27,2005-11-01
+2005-10-28,2005-05-13
+2005-10-29,2005-09-19
+2005-10-30,2005-03-09
+2005-10-31,2005-03-09
+2005-11-01,2005-11-01
+2005-11-02,2005-12-07
+2005-11-03,2005-12-07
+2005-11-04,2005-11-04
+2005-11-05,2005-11-04
+2005-11-06,2005-12-07
+2005-11-07,2005-03-09
+2005-11-08,2005-11-01
+2005-11-09,2005-11-01
+2005-11-10,2005-11-01
+2005-11-11,2005-11-01
+2005-11-12,2005-11-01
+2005-11-13,2005-11-01
+2005-11-14,2005-11-01
+2005-11-15,2005-11-01
+2005-11-16,2005-12-07
+2005-11-17,2005-12-07
+2005-11-18,2005-11-04
+2005-11-19,2005-11-04
+2005-11-20,2005-12-07
+2005-11-21,2005-12-07
+2005-11-22,2005-11-01
+2005-11-23,2005-12-07
+2005-11-24,2005-12-07
+2005-11-25,2005-11-04
+2005-11-26,2005-11-04
+2005-11-27,2005-12-07
+2005-11-28,2005-12-07
+2005-11-29,2005-12-07
+2005-11-30,2005-12-07
+2005-12-01,2005-12-07
+2005-12-02,2005-11-04
+2005-12-03,2005-11-01
+2005-12-04,2005-12-07
+2005-12-05,2005-12-07
+2005-12-06,2005-12-07
+2005-12-07,2005-12-07
+2005-12-08,2005-12-07
+2005-12-09,2005-11-01
+2005-12-10,2005-11-01
+2005-12-11,2005-12-07
+2005-12-12,2005-12-07
+2005-12-13,2005-12-07
+2005-12-14,2005-12-07
+2005-12-15,2005-12-07
+2005-12-16,2005-11-04
+2005-12-17,2005-11-01
+2005-12-18,2005-12-07
+2005-12-19,2005-12-07
+2005-12-20,2005-12-07
+2005-12-21,2005-12-07
+2005-12-22,2005-12-07
+2005-12-23,2005-12-07
+2005-12-24,2005-11-04
+2005-12-25,2005-11-04
+2005-12-26,2005-11-04
+2005-12-27,2005-12-07
+2005-12-28,2005-12-07
+2005-12-29,2005-12-07
+2005-12-30,2005-11-04
+2005-12-31,2005-11-04
diff --git a/src/calliope/example_models/urban_scale/model.yaml b/src/calliope/example_models/urban_scale/model.yaml
index e282bad6..e56c1302 100644
--- a/src/calliope/example_models/urban_scale/model.yaml
+++ b/src/calliope/example_models/urban_scale/model.yaml
@@ -13,11 +13,11 @@ config:
     calliope_version: 0.7.0
     # Time series data path - can either be a path relative to this file, or an absolute path
     time_subset: ["2005-07-01", "2005-07-02"] # Subset of timesteps
-    add_math: ["additional_math.yaml"]
 
   build:
     mode: plan # Choices: plan, operate
     ensure_feasibility: true # Switching on unmet demand
+    add_math: ["additional_math.yaml"]
 
   solve:
     solver: cbc
diff --git a/src/calliope/io.py b/src/calliope/io.py
index 4b5b4f80..205ffe7f 100644
--- a/src/calliope/io.py
+++ b/src/calliope/io.py
@@ -3,6 +3,7 @@
 """Functions to read and save model results."""
 
 import importlib.resources
+from copy import deepcopy
 from pathlib import Path
 
 # We import netCDF4 before xarray to mitigate a numpy warning:
@@ -124,16 +125,13 @@ def _deserialise(attrs: dict) -> None:
         attrs[attr] = set(attrs[attr])
 
 
-def save_netcdf(model_data, path, model=None):
+def save_netcdf(model_data, path, **kwargs):
     """Save the model to a netCDF file."""
-    original_model_data_attrs = model_data.attrs
-    model_data_attrs = original_model_data_attrs.copy()
+    original_model_data_attrs = deepcopy(model_data.attrs)
+    for key, value in kwargs.items():
+        model_data.attrs[key] = value
 
-    if model is not None and hasattr(model, "_model_def_dict"):
-        # Attach initial model definition to _model_data
-        model_data_attrs["_model_def_dict"] = model._model_def_dict.to_yaml()
-
-    _serialise(model_data_attrs)
+    _serialise(model_data.attrs)
     for var in model_data.data_vars.values():
         _serialise(var.attrs)
 
@@ -147,7 +145,6 @@ def save_netcdf(model_data, path, model=None):
     }
 
     try:
-        model_data.attrs = model_data_attrs
         model_data.to_netcdf(path, format="netCDF4", encoding=encoding)
         model_data.close()  # Force-close NetCDF file after writing
     finally:  # Revert model_data.attrs back
diff --git a/src/calliope/math/operate.yaml b/src/calliope/math/operate.yaml
index 33d3bb97..f380606f 100644
--- a/src/calliope/math/operate.yaml
+++ b/src/calliope/math/operate.yaml
@@ -1,3 +1,6 @@
+import:
+  - plan.yaml
+
 constraints:
   flow_capacity_per_storage_capacity_min.active: false
   flow_capacity_per_storage_capacity_max.active: false
diff --git a/src/calliope/math/base.yaml b/src/calliope/math/plan.yaml
similarity index 99%
rename from src/calliope/math/base.yaml
rename to src/calliope/math/plan.yaml
index 6ccc668e..0ee9f7db 100644
--- a/src/calliope/math/base.yaml
+++ b/src/calliope/math/plan.yaml
@@ -196,10 +196,10 @@ constraints:
             (
               (timesteps=get_val_at_index(timesteps=0) AND cyclic_storage=True)
               OR NOT timesteps=get_val_at_index(timesteps=0)
-            ) AND NOT lookup_cluster_first_timestep=True
+            ) AND NOT cluster_first_timestep=True
           expression: (1 - storage_loss) ** roll(timestep_resolution, timesteps=1) * roll(storage, timesteps=1)
         - where: >-
-            lookup_cluster_first_timestep=True AND NOT
+            cluster_first_timestep=True AND NOT
             (timesteps=get_val_at_index(timesteps=0) AND NOT cyclic_storage=True)
           expression: >-
             (1 - storage_loss) **
diff --git a/src/calliope/math/spores.yaml b/src/calliope/math/spores.yaml
index 74353127..28650c7d 100644
--- a/src/calliope/math/spores.yaml
+++ b/src/calliope/math/spores.yaml
@@ -1,3 +1,6 @@
+import:
+  - plan.yaml
+
 constraints:
   cost_sum_max:
     equations:
diff --git a/src/calliope/model.py b/src/calliope/model.py
index 3b1ea3db..42f67875 100644
--- a/src/calliope/model.py
+++ b/src/calliope/model.py
@@ -12,17 +12,14 @@
 import xarray as xr
 
 import calliope
-from calliope import backend, exceptions, io
-from calliope._version import __version__
+from calliope import backend, exceptions, io, preprocess
 from calliope.attrdict import AttrDict
 from calliope.postprocess import postprocess as postprocess_results
-from calliope.preprocess import load
 from calliope.preprocess.data_tables import DataTable
 from calliope.preprocess.model_data import ModelDataFactory
 from calliope.util.logging import log_time
 from calliope.util.schema import (
     CONFIG_SCHEMA,
-    MATH_SCHEMA,
     MODEL_SCHEMA,
     extract_from_schema,
     update_then_validate_config,
@@ -45,7 +42,8 @@ def read_netcdf(path):
 class Model:
     """A Calliope Model."""
 
-    _TS_OFFSET = pd.Timedelta(nanoseconds=1)
+    _TS_OFFSET = pd.Timedelta(1, unit="nanoseconds")
+    ATTRS_SAVED = ("_def_path", "applied_math")
 
     def __init__(
         self,
@@ -78,10 +76,9 @@ def __init__(
         self._timings: dict = {}
         self.config: AttrDict
         self.defaults: AttrDict
-        self.math: AttrDict
-        self._model_def_path: Path | None
+        self.applied_math: preprocess.CalliopeMath
+        self._def_path: str | None = None
         self.backend: BackendModel
-        self.math_documentation = backend.MathDocumentation()
         self._is_built: bool = False
         self._is_solved: bool = False
 
@@ -93,11 +90,16 @@ def __init__(
         if isinstance(model_definition, xr.Dataset):
             self._init_from_model_data(model_definition)
         else:
-            (model_def, self._model_def_path, applied_overrides) = (
-                load.load_model_definition(
-                    model_definition, scenario, override_dict, **kwargs
-                )
+            if isinstance(model_definition, dict):
+                model_def_dict = AttrDict(model_definition)
+            else:
+                self._def_path = str(model_definition)
+                model_def_dict = AttrDict.from_yaml(model_definition)
+
+            (model_def, applied_overrides) = preprocess.load_scenario_overrides(
+                model_def_dict, scenario, override_dict, **kwargs
             )
+
             self._init_from_model_def_dict(
                 model_def, applied_overrides, scenario, data_table_dfs
             )
@@ -111,8 +113,6 @@ def __init__(
                 f"but you are running {version_init}. Proceed with caution!"
             )
 
-        self.math_documentation.inputs = self._model_data
-
     @property
     def name(self):
         """Get the model name."""
@@ -156,7 +156,6 @@ def _init_from_model_def_dict(
         # First pass to check top-level keys are all good
         validate_dict(model_definition, CONFIG_SCHEMA, "Model definition")
 
-        self._model_def_dict = model_definition
         log_time(
             LOGGER,
             self._timings,
@@ -167,18 +166,16 @@ def _init_from_model_def_dict(
         model_config.union(model_definition.pop("config"), allow_override=True)
 
         init_config = update_then_validate_config("init", model_config)
-        # We won't store `init` in `self.config`, so we pop it out now.
-        model_config.pop("init")
 
         if init_config["time_cluster"] is not None:
             init_config["time_cluster"] = relative_path(
-                self._model_def_path, init_config["time_cluster"]
+                self._def_path, init_config["time_cluster"]
             )
 
         param_metadata = {"default": extract_from_schema(MODEL_SCHEMA, "default")}
         attributes = {
             "calliope_version_defined": init_config["calliope_version"],
-            "calliope_version_initialised": __version__,
+            "calliope_version_initialised": calliope.__version__,
             "applied_overrides": applied_overrides,
             "scenario": scenario,
             "defaults": param_metadata["default"],
@@ -186,11 +183,7 @@ def _init_from_model_def_dict(
 
         data_tables = [
             DataTable(
-                init_config,
-                source_name,
-                source_dict,
-                data_table_dfs,
-                self._model_def_path,
+                init_config, source_name, source_dict, data_table_dfs, self._def_path
             )
             for source_name, source_dict in model_definition.pop(
                 "data_tables", {}
@@ -213,9 +206,6 @@ def _init_from_model_def_dict(
 
         self._add_observed_dict("config", model_config)
 
-        math = self._add_math(init_config["add_math"])
-        self._add_observed_dict("math", math)
-
         self._model_data.attrs["name"] = init_config["name"]
         log_time(
             LOGGER,
@@ -233,11 +223,12 @@ def _init_from_model_data(self, model_data: xr.Dataset) -> None:
             model_data (xr.Dataset):
                 Model dataset with input parameters as arrays and configuration stored in the dataset attributes dictionary.
         """
-        if "_model_def_dict" in model_data.attrs:
-            self._model_def_dict = AttrDict.from_yaml_string(
-                model_data.attrs["_model_def_dict"]
+        if "_def_path" in model_data.attrs:
+            self._def_path = model_data.attrs.pop("_def_path")
+        if "applied_math" in model_data.attrs:
+            self.applied_math = preprocess.CalliopeMath.from_dict(
+                model_data.attrs.pop("applied_math")
             )
-            del model_data.attrs["_model_def_dict"]
 
         self._model_data = model_data
         self._add_model_data_methods()
@@ -260,7 +251,6 @@ def _add_model_data_methods(self):
 
         """
         self._add_observed_dict("config")
-        self._add_observed_dict("math")
 
     def _add_observed_dict(self, name: str, dict_to_add: dict | None = None) -> None:
         """Add the same dictionary as property of model object and an attribute of the model xarray dataset.
@@ -294,52 +284,18 @@ def _add_observed_dict(self, name: str, dict_to_add: dict | None = None) -> None
         self._model_data.attrs[name] = dict_to_add
         setattr(self, name, dict_to_add)
 
-    def _add_math(self, add_math: list) -> AttrDict:
-        """Load the base math and optionally override with additional math from a list of references to math files.
-
-        Args:
-            add_math (list):
-                List of references to files containing mathematical formulations that will be merged with the base formulation.
-
-        Raises:
-            exceptions.ModelError:
-                Referenced pre-defined math files or user-defined math files must exist.
-
-        Returns:
-            AttrDict: Dictionary of math (constraints, variables, objectives, and global expressions).
-        """
-        math_dir = Path(calliope.__file__).parent / "math"
-        base_math = AttrDict.from_yaml(math_dir / "base.yaml")
-
-        file_errors = []
-
-        for filename in add_math:
-            if not f"{filename}".endswith((".yaml", ".yml")):
-                yaml_filepath = math_dir / f"{filename}.yaml"
-            else:
-                yaml_filepath = relative_path(self._model_def_path, filename)
-
-            if not yaml_filepath.is_file():
-                file_errors.append(filename)
-                continue
-            else:
-                override_dict = AttrDict.from_yaml(yaml_filepath)
-
-            base_math.union(override_dict, allow_override=True)
-        if file_errors:
-            raise exceptions.ModelError(
-                f"Attempted to load additional math that does not exist: {file_errors}"
-            )
-        self._model_data.attrs["applied_additional_math"] = add_math
-        return base_math
-
-    def build(self, force: bool = False, **kwargs) -> None:
+    def build(
+        self, force: bool = False, add_math_dict: dict | None = None, **kwargs
+    ) -> None:
         """Build description of the optimisation problem in the chosen backend interface.
 
         Args:
             force (bool, optional):
                 If ``force`` is True, any existing results will be overwritten.
                 Defaults to False.
+            add_math_dict (dict | None, optional):
+                Additional math to apply on top of the YAML base / additional math files.
+                Content of this dictionary will override any matching key:value pairs in the loaded math files.
             **kwargs: build configuration overrides.
         """
         if self._is_built and not force:
@@ -355,7 +311,8 @@ def build(self, force: bool = False, **kwargs) -> None:
         )
 
         backend_config = {**self.config["build"], **kwargs}
-        if backend_config["mode"] == "operate":
+        mode = backend_config["mode"]
+        if mode == "operate":
             if not self._model_data.attrs["allow_operate_mode"]:
                 raise exceptions.ModelError(
                     "Unable to run this model in operate (i.e. dispatch) mode, probably because "
@@ -367,11 +324,20 @@ def build(self, force: bool = False, **kwargs) -> None:
             )
         else:
             backend_input = self._model_data
+
+        init_math_list = [] if backend_config.get("ignore_mode_math") else [mode]
+        end_math_list = [] if add_math_dict is None else [add_math_dict]
+        full_math_list = init_math_list + backend_config["add_math"] + end_math_list
+        LOGGER.debug(f"Math preprocessing | Loading math: {full_math_list}")
+        model_math = preprocess.CalliopeMath(full_math_list, self._def_path)
+
         backend_name = backend_config.pop("backend")
         self.backend = backend.get_model_backend(
-            backend_name, backend_input, **backend_config
+            backend_name, backend_input, model_math, **backend_config
         )
-        self.backend.add_all_math()
+        self.backend.add_optimisation_components()
+
+        self.applied_math = model_math
 
         self._model_data.attrs["timestamp_build_complete"] = log_time(
             LOGGER,
@@ -497,7 +463,14 @@ def run(self, force_rerun=False, **kwargs):
 
     def to_netcdf(self, path):
         """Save complete model data (inputs and, if available, results) to a NetCDF file at the given `path`."""
-        io.save_netcdf(self._model_data, path, model=self)
+        saved_attrs = {}
+        for attr in set(self.ATTRS_SAVED) & set(self.__dict__.keys()):
+            if not isinstance(getattr(self, attr), str | list | None):
+                saved_attrs[attr] = dict(getattr(self, attr))
+            else:
+                saved_attrs[attr] = getattr(self, attr)
+
+        io.save_netcdf(self._model_data, path, **saved_attrs)
 
     def to_csv(
         self, path: str | Path, dropna: bool = True, allow_overwrite: bool = False
@@ -533,60 +506,6 @@ def info(self) -> str:
         )
         return "\n".join(info_strings)
 
-    def validate_math_strings(self, math_dict: dict) -> None:
-        """Validate that `expression` and `where` strings of a dictionary containing string mathematical formulations can be successfully parsed.
-
-        This function can be used to test user-defined math before attempting to build the optimisation problem.
-
-        NOTE: strings are not checked for evaluation validity. Evaluation issues will be raised only on calling `Model.build()`.
-
-        Args:
-            math_dict (dict): Math formulation dictionary to validate. Top level keys must be one or more of ["variables", "global_expressions", "constraints", "objectives"], e.g.:
-                ```python
-                {
-                    "constraints": {
-                        "my_constraint_name":
-                            {
-                                "foreach": ["nodes"],
-                                "where": "base_tech=supply",
-                                "equations": [{"expression": "sum(flow_cap, over=techs) >= 10"}]
-                            }
-
-                        }
-                }
-                ```
-        Returns:
-            If all components of the dictionary are parsed successfully, this function will log a success message to the INFO logging level and return None.
-            Otherwise, a calliope.ModelError will be raised with parsing issues listed.
-        """
-        validate_dict(math_dict, MATH_SCHEMA, "math")
-        valid_component_names = [
-            *self.math["variables"].keys(),
-            *self.math["global_expressions"].keys(),
-            *math_dict.get("variables", {}).keys(),
-            *math_dict.get("global_expressions", {}).keys(),
-            *self.inputs.data_vars.keys(),
-            *self.inputs.attrs["defaults"].keys(),
-        ]
-        collected_errors: dict = dict()
-        for component_group, component_dicts in math_dict.items():
-            for name, component_dict in component_dicts.items():
-                parsed = backend.ParsedBackendComponent(
-                    component_group, name, component_dict
-                )
-                parsed.parse_top_level_where(errors="ignore")
-                parsed.parse_equations(set(valid_component_names), errors="ignore")
-                if not parsed._is_valid:
-                    collected_errors[f"{component_group}:{name}"] = parsed._errors
-
-        if collected_errors:
-            exceptions.print_warnings_and_raise_errors(
-                during="math string parsing (marker indicates where parsing stopped, which might not be the root cause of the issue; sorry...)",
-                errors=collected_errors,
-            )
-
-        LOGGER.info("Model: validated math strings")
-
     def _prepare_operate_mode_inputs(
         self, start_window_idx: int = 0, **config_kwargs
     ) -> xr.Dataset:
diff --git a/src/calliope/postprocess/math_documentation.py b/src/calliope/postprocess/math_documentation.py
new file mode 100644
index 00000000..ebfb3193
--- /dev/null
+++ b/src/calliope/postprocess/math_documentation.py
@@ -0,0 +1,101 @@
+"""Post-processing functions to create math documentation."""
+
+import logging
+import typing
+from pathlib import Path
+from typing import Literal, overload
+
+from calliope.backend import ALLOWED_MATH_FILE_FORMATS, LatexBackendModel
+from calliope.model import Model
+
+LOGGER = logging.getLogger(__name__)
+
+
+class MathDocumentation:
+    """For math documentation."""
+
+    def __init__(
+        self, model: Model, include: Literal["all", "valid"] = "all", **kwargs
+    ) -> None:
+        """Math documentation builder/writer.
+
+        Backend is always built by default.
+
+        Args:
+            model (Model): initialised Callipe model instance.
+            include (Literal["all", "valid"], optional):
+                Either include all possible math equations ("all") or only those for
+                which at least one "where" case is valid ("valid"). Defaults to "all".
+            **kwargs: kwargs for the LaTeX backend.
+        """
+        self.name: str = model.name + " math"
+        self.backend: LatexBackendModel = LatexBackendModel(
+            model._model_data, model.applied_math, include, **kwargs
+        )
+        self.backend.add_optimisation_components()
+
+    @property
+    def math(self):
+        """Direct access to backend math."""
+        return self.backend.math
+
+    # Expecting string if not giving filename.
+    @overload
+    def write(
+        self,
+        filename: Literal[None] = None,
+        mkdocs_features: bool = False,
+        format: ALLOWED_MATH_FILE_FORMATS | None = None,
+    ) -> str: ...
+
+    # Expecting None (and format arg is not needed) if giving filename.
+    @overload
+    def write(self, filename: str | Path, mkdocs_features: bool = False) -> None: ...
+
+    def write(
+        self,
+        filename: str | Path | None = None,
+        mkdocs_features: bool = False,
+        format: ALLOWED_MATH_FILE_FORMATS | None = None,
+    ) -> str | None:
+        """Write model documentation.
+
+        Args:
+            filename (str | Path | None, optional):
+                If given, will write the built mathematical formulation to a file with
+                the given extension as the file format. Defaults to None.
+            mkdocs_features (bool, optional):
+                If True and Markdown docs are being generated, then:
+                - the equations will be on a tab and the original YAML math definition will be on another tab;
+                - the equation cross-references will be given in a drop-down list.
+                Defaults to False.
+            format (ALLOWED_MATH_FILE_FORMATS | None, optional):
+                Not required if filename is given (as the format will be automatically inferred).
+                Required if expecting a string return from calling this function. The LaTeX math will be embedded in a document of the given format (tex=LaTeX, rst=reStructuredText, md=Markdown).
+                Defaults to None.
+
+        Raises:
+            ValueError: The file format (inferred automatically from `filename` or given by `format`) must be one of ["tex", "rst", "md"].
+
+        Returns:
+            str | None:
+                If `filename` is None, the built mathematical formulation documentation will be returned as a string.
+        """
+        if format is None and filename is not None:
+            format = Path(filename).suffix.removeprefix(".")  # type: ignore
+            LOGGER.info(
+                f"Inferring math documentation format from filename as `{format}`."
+            )
+
+        allowed_formats = typing.get_args(ALLOWED_MATH_FILE_FORMATS)
+        if format is None or format not in allowed_formats:
+            raise ValueError(
+                f"Math documentation format must be one of {allowed_formats}, received `{format}`"
+            )
+        populated_doc = self.backend.generate_math_doc(format, mkdocs_features)
+
+        if filename is None:
+            return populated_doc
+        else:
+            Path(filename).write_text(populated_doc)
+            return None
diff --git a/src/calliope/preprocess/__init__.py b/src/calliope/preprocess/__init__.py
index 0ba2ad52..2b9584be 100644
--- a/src/calliope/preprocess/__init__.py
+++ b/src/calliope/preprocess/__init__.py
@@ -1 +1,6 @@
 """Preprocessing module."""
+
+from calliope.preprocess.data_tables import DataTable
+from calliope.preprocess.model_data import ModelDataFactory
+from calliope.preprocess.model_math import CalliopeMath
+from calliope.preprocess.scenarios import load_scenario_overrides
diff --git a/src/calliope/preprocess/model_math.py b/src/calliope/preprocess/model_math.py
new file mode 100644
index 00000000..a05a6a12
--- /dev/null
+++ b/src/calliope/preprocess/model_math.py
@@ -0,0 +1,176 @@
+"""Calliope math handling with interfaces for pre-defined and user-defined files."""
+
+import importlib.resources
+import logging
+import typing
+from copy import deepcopy
+from pathlib import Path
+
+from calliope.attrdict import AttrDict
+from calliope.exceptions import ModelError
+from calliope.util.schema import MATH_SCHEMA, validate_dict
+from calliope.util.tools import relative_path
+
+LOGGER = logging.getLogger(__name__)
+ORDERED_COMPONENTS_T = typing.Literal[
+    "variables",
+    "global_expressions",
+    "constraints",
+    "piecewise_constraints",
+    "objectives",
+]
+
+
+class CalliopeMath:
+    """Calliope math handling."""
+
+    ATTRS_TO_SAVE = ("history", "data")
+    ATTRS_TO_LOAD = ("history",)
+
+    def __init__(
+        self, math_to_add: list[str | dict], model_def_path: str | Path | None = None
+    ):
+        """Calliope YAML math handler.
+
+        Args:
+            math_to_add (list[str | dict]):
+                List of Calliope math to load.
+                If a string, it can be a reference to pre-/user-defined math files.
+                If a dictionary, it is equivalent in structure to a YAML math file.
+            model_def_path (str | Path | None, optional): Model definition path, needed when using relative paths. Defaults to None.
+        """
+        self.history: list[str] = []
+        self.data: AttrDict = AttrDict(
+            {name: {} for name in typing.get_args(ORDERED_COMPONENTS_T)}
+        )
+
+        for math in math_to_add:
+            if isinstance(math, dict):
+                self.add(AttrDict(math))
+            else:
+                self._init_from_string(math, model_def_path)
+
+    def __eq__(self, other):
+        """Compare between two model math instantiations."""
+        if not isinstance(other, CalliopeMath):
+            return NotImplemented
+        return self.history == other.history and self.data == other.data
+
+    def __iter__(self):
+        """Enable dictionary conversion."""
+        for key in self.ATTRS_TO_SAVE:
+            yield key, deepcopy(getattr(self, key))
+
+    def __repr__(self) -> str:
+        """Custom string representation of class."""
+        return f"""Calliope math definition dictionary with:
+    {len(self.data["variables"])} decision variable(s)
+    {len(self.data["global_expressions"])} global expression(s)
+    {len(self.data["constraints"])} constraint(s)
+    {len(self.data["piecewise_constraints"])} piecewise constraint(s)
+    {len(self.data["objectives"])} objective(s)
+        """
+
+    def add(self, math: AttrDict):
+        """Add math into the model.
+
+        Args:
+            math (AttrDict): Valid math dictionary.
+        """
+        self.data.union(math, allow_override=True)
+
+    @classmethod
+    def from_dict(cls, math_dict: dict) -> "CalliopeMath":
+        """Load a CalliopeMath object from a dictionary representation, recuperating relevant attributes.
+
+        Args:
+            math_dict (dict): Dictionary representation of a CalliopeMath object.
+
+        Returns:
+            CalliopeMath: Loaded from supplied dictionary representation.
+        """
+        new_self = cls([math_dict["data"]])
+        for attr in cls.ATTRS_TO_LOAD:
+            setattr(new_self, attr, math_dict[attr])
+        return new_self
+
+    def in_history(self, math_name: str) -> bool:
+        """Evaluate if math has already been applied.
+
+        Args:
+            math_name (str): Math file to check.
+
+        Returns:
+            bool: `True` if found in history. `False` otherwise.
+        """
+        return math_name in self.history
+
+    def validate(self) -> None:
+        """Test current math and optional external math against the MATH schema."""
+        validate_dict(self.data, MATH_SCHEMA, "math")
+        LOGGER.info("Math preprocessing | validated math against schema.")
+
+    def _add_pre_defined_file(self, filename: str) -> None:
+        """Add pre-defined Calliope math.
+
+        Args:
+            filename (str): name of Calliope internal math (no suffix).
+
+        Raises:
+            ModelError: If math has already been applied.
+        """
+        if self.in_history(filename):
+            raise ModelError(
+                f"Math preprocessing | Overwriting with previously applied pre-defined math: '{filename}'."
+            )
+        with importlib.resources.as_file(
+            importlib.resources.files("calliope") / "math"
+        ) as f:
+            self._add_file(f / f"{filename}.yaml", filename)
+
+    def _add_user_defined_file(
+        self, relative_filepath: str | Path, model_def_path: str | Path | None
+    ) -> None:
+        """Add user-defined Calliope math, relative to the model definition path.
+
+        Args:
+            relative_filepath (str | Path): Path to user math, relative to model definition.
+            model_def_path (str | Path): Model definition path.
+
+        Raises:
+            ModelError: If file has already been applied.
+        """
+        math_name = str(relative_filepath)
+        if self.in_history(math_name):
+            raise ModelError(
+                f"Math preprocessing | Overwriting with previously applied user-defined math: '{relative_filepath}'."
+            )
+        self._add_file(relative_path(model_def_path, relative_filepath), math_name)
+
+    def _init_from_string(
+        self, math_to_add: str, model_def_path: str | Path | None = None
+    ):
+        """Load math definition from a list of files.
+
+        Args:
+            math_to_add (str): Calliope math file to load. Suffix implies user-math.
+            model_def_path (str | Path | None, optional): Model definition path. Defaults to None.
+
+        Raises:
+            ModelError: User-math requested without providing `model_def_path`.
+        """
+        if not math_to_add.endswith((".yaml", ".yml")):
+            self._add_pre_defined_file(math_to_add)
+        else:
+            self._add_user_defined_file(math_to_add, model_def_path)
+
+    def _add_file(self, yaml_filepath: Path, name: str) -> None:
+        try:
+            math = AttrDict.from_yaml(yaml_filepath, allow_override=True)
+        except FileNotFoundError:
+            raise ModelError(
+                f"Math preprocessing | File does not exist: {yaml_filepath}"
+            )
+        self.add(math)
+        self.history.append(name)
+        LOGGER.info(f"Math preprocessing | added file '{name}'.")
diff --git a/src/calliope/preprocess/load.py b/src/calliope/preprocess/scenarios.py
similarity index 55%
rename from src/calliope/preprocess/load.py
rename to src/calliope/preprocess/scenarios.py
index 6c1df045..473544fb 100644
--- a/src/calliope/preprocess/load.py
+++ b/src/calliope/preprocess/scenarios.py
@@ -3,7 +3,6 @@
 """Preprocessing of base model definition and overrides/scenarios into a unified dictionary."""
 
 import logging
-from pathlib import Path
 
 from calliope import exceptions
 from calliope.attrdict import AttrDict
@@ -12,96 +11,35 @@
 LOGGER = logging.getLogger(__name__)
 
 
-def load_model_definition(
-    model_definition: str | Path | dict,
+def load_scenario_overrides(
+    model_definition: dict,
     scenario: str | None = None,
     override_dict: dict | None = None,
     **kwargs,
-) -> tuple[AttrDict, Path | None, str]:
-    """Load model definition from file / dictionary and apply user-defined overrides.
+) -> tuple[AttrDict, str]:
+    """Apply user-defined overrides to the model definition.
 
     Args:
-        model_definition (str | Path | dict):
-            If string or pathlib.Path, path to YAML file with model configuration.
-            If dictionary, equivalent to loading the model configuration YAML from file.
-        scenario (str | None, optional):
-            If not None, name of scenario to apply.
-            Can either be a named scenario, or a comma-separated list of individual overrides to be combined ad-hoc,
-            e.g. 'my_scenario_name' or 'override1,override2'.
+        model_definition (dict):
+            Model definition dictionary.
+        scenario (str | None, optional): Scenario(s) to apply, comma separated.
+            e.g.: 'my_scenario_name' or 'override1,override2'.
             Defaults to None.
         override_dict (dict | None, optional):
-            If not None, dictionary of overrides to apply.
-            These will be applied _after_ `scenario` overrides.
+            Overrides to apply _after_ `scenario` overrides.
             Defaults to None.
-        **kwargs: initialisation overrides.
+        **kwargs:
+            initialisation overrides.
 
     Returns:
-        tuple[AttrDict, Path | None, str]:
+        tuple[AttrDict, str]:
             1. Model definition with overrides applied.
-            2. Path to model definition YAML if input `model_definiton` was pathlike, otherwise None.
-            3. Expansion of scenarios (which are references to model overrides) into a list of named override(s) that have been applied.
+            2. Expansion of scenarios (which are references to model overrides) into a list of named override(s) that have been applied.
     """
-    if not isinstance(model_definition, dict):
-        model_def_path = Path(model_definition)
-        model_def_dict = AttrDict.from_yaml(model_def_path)
-    else:
-        model_def_dict = AttrDict(model_definition)
-        model_def_path = None
-
-    model_def_with_overrides, applied_overrides = _apply_overrides(
-        model_def_dict, scenario=scenario, override_dict=override_dict
-    )
-    model_def_with_overrides.union(
-        AttrDict({"config.init": kwargs}), allow_override=True
-    )
-
-    return (model_def_with_overrides, model_def_path, ";".join(applied_overrides))
-
-
-def _combine_overrides(overrides: AttrDict, scenario_overrides: list):
-    combined_override_dict = AttrDict()
-    for override in scenario_overrides:
-        try:
-            yaml_string = overrides[override].to_yaml()
-            override_with_imports = AttrDict.from_yaml_string(yaml_string)
-        except KeyError:
-            raise exceptions.ModelError(f"Override `{override}` is not defined.")
-        try:
-            combined_override_dict.union(override_with_imports, allow_override=False)
-        except KeyError as e:
-            raise exceptions.ModelError(
-                f"{str(e)[1:-1]}. Already specified but defined again in override `{override}`."
-            )
-
-    return combined_override_dict
-
-
-def _apply_overrides(
-    model_def: AttrDict,
-    scenario: str | None = None,
-    override_dict: str | dict | None = None,
-) -> tuple[AttrDict, list[str]]:
-    """Generate processed Model configuration, applying any scenario overrides.
-
-    Args:
-        model_def (calliope.Attrdict): Loaded model definition as an attribute dictionary.
-        scenario (str | None, optional):
-            If not None, name of scenario to apply.
-            Can either be a named scenario, or a comma-separated list of individual overrides to be combined ad-hoc,
-            e.g. 'my_scenario_name' or 'override1,override2'.
-            Defaults to None.
-        override_dict (str | dict | None, optional):
-            If not None, dictionary of overrides to apply.
-            These will be applied _after_ `scenario` overrides.
-            Defaults to None.
+    model_def_dict = AttrDict(model_definition)
 
-    Returns:
-        tuple[AttrDict, list[str]]:
-            1. Model definition dictionary with overrides applied from `scenario` and `override_dict`.
-            1. Expansion of scenarios (which are references to model overrides) into a list of named override(s) that have been applied.
-    """
     # The input files are allowed to override other model defaults
-    model_def_copy = model_def.copy()
+    model_def_with_overrides = model_def_dict.copy()
 
     # First pass of applying override dict before applying scenarios,
     # so that can override scenario definitions by override_dict
@@ -110,43 +48,69 @@ def _apply_overrides(
 
     if isinstance(override_dict, dict):
         override_dict = AttrDict(override_dict)
-        model_def_copy.union(override_dict, allow_override=True, allow_replacement=True)
+        model_def_with_overrides.union(
+            override_dict, allow_override=True, allow_replacement=True
+        )
 
-    overrides = model_def_copy.pop("overrides", {})
-    scenarios = model_def_copy.pop("scenarios", {})
+    overrides = model_def_with_overrides.pop("overrides", {})
+    scenarios = model_def_with_overrides.pop("scenarios", {})
 
     if scenario is not None:
-        scenario_overrides = _load_overrides_from_scenario(
-            model_def_copy, scenario, overrides, scenarios
+        applied_overrides = _load_overrides_from_scenario(
+            model_def_with_overrides, scenario, overrides, scenarios
         )
         LOGGER.info(
-            f"(scenarios, {scenario} ) | Applying the following overrides: {scenario_overrides}."
+            f"(scenarios, {scenario} ) | Applying the following overrides: {applied_overrides}."
         )
-        overrides_from_scenario = _combine_overrides(overrides, scenario_overrides)
+        overrides_from_scenario = _combine_overrides(overrides, applied_overrides)
 
-        model_def_copy.union(
+        model_def_with_overrides.union(
             overrides_from_scenario, allow_override=True, allow_replacement=True
         )
     else:
-        scenario_overrides = []
+        applied_overrides = []
 
     # Second pass of applying override dict after applying scenarios,
     # so that scenario-based overrides are overridden by override_dict!
     if override_dict is not None:
-        model_def_copy.union(override_dict, allow_override=True, allow_replacement=True)
-    if "locations" in model_def_copy.keys():
+        model_def_with_overrides.union(
+            override_dict, allow_override=True, allow_replacement=True
+        )
+    if "locations" in model_def_with_overrides.keys():
         # TODO: remove in v0.7.1
         exceptions.warn(
             "`locations` has been renamed to `nodes` and will stop working "
             "in v0.7.1. Please update your model configuration accordingly.",
             FutureWarning,
         )
-        model_def_copy["nodes"] = model_def_copy["locations"]
-        del model_def_copy["locations"]
+        model_def_with_overrides["nodes"] = model_def_with_overrides["locations"]
+        del model_def_with_overrides["locations"]
 
-    _log_overrides(model_def, model_def_copy)
+    _log_overrides(model_def_dict, model_def_with_overrides)
 
-    return model_def_copy, scenario_overrides
+    model_def_with_overrides.union(
+        AttrDict({"config.init": kwargs}), allow_override=True
+    )
+
+    return (model_def_with_overrides, ";".join(applied_overrides))
+
+
+def _combine_overrides(overrides: AttrDict, scenario_overrides: list):
+    combined_override_dict = AttrDict()
+    for override in scenario_overrides:
+        try:
+            yaml_string = overrides[override].to_yaml()
+            override_with_imports = AttrDict.from_yaml_string(yaml_string)
+        except KeyError:
+            raise exceptions.ModelError(f"Override `{override}` is not defined.")
+        try:
+            combined_override_dict.union(override_with_imports, allow_override=False)
+        except KeyError as e:
+            raise exceptions.ModelError(
+                f"{str(e)[1:-1]}. Already specified but defined again in override `{override}`."
+            )
+
+    return combined_override_dict
 
 
 def _load_overrides_from_scenario(
diff --git a/src/calliope/preprocess/time.py b/src/calliope/preprocess/time.py
index ebb16ee7..6b1cc97b 100644
--- a/src/calliope/preprocess/time.py
+++ b/src/calliope/preprocess/time.py
@@ -262,11 +262,14 @@ def _lookup_clusters(dataset: xr.Dataset, grouper: pd.Series) -> xr.Dataset:
     1. the first and last timestep of the cluster,
     2. the last timestep of the cluster corresponding to a date in the original timeseries
     """
-    dataset["lookup_cluster_first_timestep"] = dataset.timesteps.isin(
+    dataset["cluster_first_timestep"] = dataset.timesteps.isin(
         dataset.timesteps.groupby("timesteps.date").first()
     )
-    dataset["lookup_cluster_last_timestep"] = dataset.timesteps.isin(
-        dataset.timesteps.groupby("timesteps.date").last()
+    dataset["lookup_cluster_last_timestep"] = (
+        dataset.timesteps.groupby("timesteps.date")
+        .last()
+        .rename({"date": "timesteps"})
+        .reindex_like(dataset.timesteps)
     )
 
     dataset["lookup_datestep_cluster"] = xr.DataArray(
diff --git a/src/calliope/util/tools.py b/src/calliope/util/tools.py
index 82194927..51920d88 100644
--- a/src/calliope/util/tools.py
+++ b/src/calliope/util/tools.py
@@ -19,13 +19,13 @@ def relative_path(base_path_file, path) -> Path:
     """
     # Check if base_path_file is a string because it might be an AttrDict
     path = Path(path)
-    if base_path_file is not None:
+    if path.is_absolute() or base_path_file is None:
+        return path
+    else:
         base_path_file = Path(base_path_file)
         if base_path_file.is_file():
             base_path_file = base_path_file.parent
-        if not path.is_absolute():
-            path = base_path_file.absolute() / path
-    return path
+        return base_path_file.absolute() / path
 
 
 def listify(var: Any) -> list:
@@ -40,7 +40,9 @@ def listify(var: Any) -> list:
     Returns:
         list: List containing `var` or elements of `var` (if input was a non-string iterable).
     """
-    if not isinstance(var, str) and hasattr(var, "__iter__"):
+    if var is None:
+        var = []
+    elif not isinstance(var, str) and hasattr(var, "__iter__"):
         var = list(var)
     else:
         var = [var]
diff --git a/tests/common/util.py b/tests/common/util.py
index 5f0fe356..8ae70da8 100644
--- a/tests/common/util.py
+++ b/tests/common/util.py
@@ -5,7 +5,8 @@
 import xarray as xr
 
 import calliope
-from calliope import backend
+import calliope.backend
+import calliope.preprocess
 
 
 def build_test_model(
@@ -80,9 +81,9 @@ def check_variable_exists(
 def build_lp(
     model: calliope.Model,
     outfile: str | Path,
-    math: dict[str, dict | list] | None = None,
+    math_data: dict[str, dict | list] | None = None,
     backend_name: Literal["pyomo"] = "pyomo",
-) -> "backend.BackendModel":
+) -> "calliope.backend.backend_model.BackendModel":
     """
     Write a barebones LP file with which to compare in tests.
     All model parameters and variables will be loaded automatically, as well as a dummy objective if one isn't provided as part of `math`.
@@ -94,41 +95,43 @@ def build_lp(
         math (dict | None, optional): All constraint/global expression/objective math to apply. Defaults to None.
         backend_name (Literal["pyomo"], optional): Backend to use to create the LP file. Defaults to "pyomo".
     """
-    backend_instance = backend.get_model_backend(backend_name, model._model_data)
-
-    for name, dict_ in model.math["variables"].items():
-        backend_instance.add_variable(name, dict_)
-    for name, dict_ in model.math["global_expressions"].items():
-        backend_instance.add_global_expression(name, dict_)
+    math = calliope.preprocess.CalliopeMath(
+        ["plan", *model.config.build.get("add_math", [])]
+    )
 
-    if isinstance(math, dict):
-        for component_group, component_math in math.items():
-            component = component_group.removesuffix("s")
+    math_to_add = calliope.AttrDict()
+    if isinstance(math_data, dict):
+        for component_group, component_math in math_data.items():
             if isinstance(component_math, dict):
-                for name, dict_ in component_math.items():
-                    getattr(backend_instance, f"add_{component}")(name, dict_)
+                math_to_add.union(calliope.AttrDict({component_group: component_math}))
             elif isinstance(component_math, list):
                 for name in component_math:
-                    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():
-        backend_instance.add_objective(
-            "dummy_obj", {"equations": [{"expression": "1 + 1"}], "sense": "minimize"}
-        )
-        backend_instance._instance.objectives["dummy_obj"][0].activate()
-    elif "objectives" in math.keys():
-        if isinstance(math["objectives"], dict):
-            objectives = list(math["objectives"].keys())
-        else:
-            objectives = math["objectives"]
-        assert len(objectives) == 1, "Can only test with one objective"
-        backend_instance._instance.objectives[objectives[0]][0].activate()
+                    math_to_add.set_key(
+                        f"{component_group}.{name}", math.data[component_group][name]
+                    )
+    if math_data is None or "objectives" not in math_to_add.keys():
+        obj = {
+            "dummy_obj": {"equations": [{"expression": "1 + 1"}], "sense": "minimize"}
+        }
+        math_to_add.union(calliope.AttrDict({"objectives": obj}))
+        obj_to_activate = "dummy_obj"
+    else:
+        obj_to_activate = list(math_to_add["objectives"].keys())[0]
+    del math.data["constraints"]
+    del math.data["objectives"]
+    math.add(math_to_add)
+
+    model.build(
+        add_math_dict=math.data,
+        ignore_mode_math=True,
+        objective=obj_to_activate,
+        add_math=[],
+        pre_validate_math_strings=False,
+    )
 
-    backend_instance.verbose_strings()
+    model.backend.verbose_strings()
 
-    backend_instance.to_lp(str(outfile))
+    model.backend.to_lp(str(outfile))
 
     # strip trailing whitespace from `outfile` after the fact,
     # so it can be reliably compared other files in future
@@ -139,4 +142,4 @@ def build_lp(
 
     # reintroduce the trailing newline since both Pyomo and file formatters love them.
     Path(outfile).write_text("\n".join(stripped_lines) + "\n")
-    return backend_instance
+    return model.backend
diff --git a/tests/conftest.py b/tests/conftest.py
index b01dc5ba..3d4694c5 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -7,6 +7,7 @@
 
 from calliope.attrdict import AttrDict
 from calliope.backend import latex_backend_model, pyomo_backend_model
+from calliope.preprocess import CalliopeMath
 from calliope.util.schema import CONFIG_SCHEMA, MODEL_SCHEMA, extract_from_schema
 
 from .common.util import build_test_model as build_model
@@ -156,6 +157,20 @@ def simple_conversion_plus():
     return m
 
 
+@pytest.fixture(scope="module")
+def dummy_model_math():
+    math = {
+        "data": {
+            "constraints": {},
+            "variables": {},
+            "global_expressions": {},
+            "objectives": {},
+        },
+        "history": [],
+    }
+    return CalliopeMath.from_dict(math)
+
+
 @pytest.fixture(scope="module")
 def dummy_model_data(config_defaults, model_defaults):
     coords = {
@@ -242,10 +257,6 @@ def dummy_model_data(config_defaults, model_defaults):
                 ["nodes", "techs"],
                 [[False, False, False, False], [False, False, False, True]],
             ),
-            "primary_carrier_out": (
-                ["carriers", "techs"],
-                [[1.0, np.nan, 1.0, np.nan], [np.nan, 1.0, np.nan, np.nan]],
-            ),
             "lookup_techs": (["techs"], ["foobar", np.nan, "foobaz", np.nan]),
             "lookup_techs_no_match": (["techs"], ["foo", np.nan, "bar", np.nan]),
             "lookup_multi_dim_nodes": (
@@ -296,13 +307,11 @@ def dummy_model_data(config_defaults, model_defaults):
     # This value is set on the parameter directly to ensure it finds its way through to the LaTex math.
     model_data.no_dims.attrs["default"] = 0
 
-    model_data.attrs["math"] = AttrDict(
-        {"constraints": {}, "variables": {}, "global_expressions": {}, "objectives": {}}
-    )
     return model_data
 
 
 def populate_backend_model(backend):
+    backend._add_all_inputs_as_parameters()
     backend.add_variable(
         "multi_dim_var",
         {
@@ -335,18 +344,20 @@ def populate_backend_model(backend):
 
 
 @pytest.fixture(scope="module")
-def dummy_pyomo_backend_model(dummy_model_data):
-    backend = pyomo_backend_model.PyomoBackendModel(dummy_model_data)
+def dummy_pyomo_backend_model(dummy_model_data, dummy_model_math):
+    backend = pyomo_backend_model.PyomoBackendModel(dummy_model_data, dummy_model_math)
     return populate_backend_model(backend)
 
 
 @pytest.fixture(scope="module")
-def dummy_latex_backend_model(dummy_model_data):
-    backend = latex_backend_model.LatexBackendModel(dummy_model_data)
+def dummy_latex_backend_model(dummy_model_data, dummy_model_math):
+    backend = latex_backend_model.LatexBackendModel(dummy_model_data, dummy_model_math)
     return populate_backend_model(backend)
 
 
 @pytest.fixture(scope="class")
-def valid_latex_backend(dummy_model_data):
-    backend = latex_backend_model.LatexBackendModel(dummy_model_data, include="valid")
+def valid_latex_backend(dummy_model_data, dummy_model_math):
+    backend = latex_backend_model.LatexBackendModel(
+        dummy_model_data, dummy_model_math, include="valid"
+    )
     return populate_backend_model(backend)
diff --git a/tests/test_backend_general.py b/tests/test_backend_general.py
index 0d634ef6..f61443f5 100644
--- a/tests/test_backend_general.py
+++ b/tests/test_backend_general.py
@@ -29,7 +29,7 @@ def built_model_cls_longnames(backend) -> calliope.Model:
 @pytest.fixture
 def built_model_func_longnames(backend) -> calliope.Model:
     m = build_model({}, "simple_supply,two_hours,investment_costs")
-    m.build(backend=backend)
+    m.build(backend=backend, pre_validate_math_strings=False)
     m.backend.verbose_strings()
     return m
 
@@ -37,7 +37,7 @@ def built_model_func_longnames(backend) -> calliope.Model:
 @pytest.fixture
 def solved_model_func(backend) -> calliope.Model:
     m = build_model({}, "simple_supply,two_hours,investment_costs")
-    m.build(backend=backend)
+    m.build(backend=backend, pre_validate_math_strings=False)
     m.solve()
     return m
 
@@ -69,7 +69,7 @@ def solved_model_cls(backend) -> calliope.Model:
 @pytest.fixture
 def built_model_func_updated_cost_flow_cap(backend, dummy_int: int) -> calliope.Model:
     m = build_model({}, "simple_supply,two_hours,investment_costs")
-    m.build(backend=backend)
+    m.build(backend=backend, pre_validate_math_strings=False)
     m.backend.verbose_strings()
     m.backend.update_parameter("cost_flow_cap", dummy_int)
     return m
diff --git a/tests/test_backend_gurobi.py b/tests/test_backend_gurobi.py
index 18b8c8eb..95821860 100755
--- a/tests/test_backend_gurobi.py
+++ b/tests/test_backend_gurobi.py
@@ -28,7 +28,7 @@ def simple_supply_gurobi(self):
     @pytest.fixture
     def simple_supply_gurobi_func(self):
         m = build_model({}, "simple_supply,two_hours,investment_costs")
-        m.build(backend="gurobi")
+        m.build(backend="gurobi", pre_validate_math_strings=False)
         m.solve()
         return m
 
diff --git a/tests/test_backend_latex_backend.py b/tests/test_backend_latex_backend.py
index e309484e..e28b0830 100644
--- a/tests/test_backend_latex_backend.py
+++ b/tests/test_backend_latex_backend.py
@@ -1,5 +1,4 @@
 import textwrap
-from pathlib import Path
 
 import pytest
 import xarray as xr
@@ -7,65 +6,7 @@
 from calliope import exceptions
 from calliope.backend import latex_backend_model
 
-from .common.util import build_test_model, check_error_or_warning
-
-
-class TestMathDocumentation:
-    @pytest.fixture(scope="class")
-    def no_build(self):
-        return build_test_model({}, "simple_supply,two_hours,investment_costs")
-
-    @pytest.fixture(scope="class")
-    def build_all(self):
-        model = build_test_model({}, "simple_supply,two_hours,investment_costs")
-        model.math_documentation.build(include="all")
-        return model
-
-    @pytest.fixture(scope="class")
-    def build_valid(self):
-        model = build_test_model({}, "simple_supply,two_hours,investment_costs")
-        model.math_documentation.build(include="valid")
-        return model
-
-    def test_write_before_build(self, no_build, tmpdir_factory):
-        filepath = tmpdir_factory.mktemp("custom_math").join("foo.tex")
-        with pytest.raises(exceptions.ModelError) as excinfo:
-            no_build.math_documentation.write(filepath)
-        assert check_error_or_warning(
-            excinfo, "Build the documentation (`build`) before trying to write it"
-        )
-
-    @pytest.mark.parametrize(
-        ("format", "startswith"),
-        [
-            ("tex", "\n\\documentclass{article}"),
-            ("rst", "\nObjective"),
-            ("md", "\n## Objective"),
-        ],
-    )
-    @pytest.mark.parametrize("include", ["build_all", "build_valid"])
-    def test_string_return(self, request, format, startswith, include):
-        model = request.getfixturevalue(include)
-        string_math = model.math_documentation.write(format=format)
-        assert string_math.startswith(startswith)
-
-    def test_to_file(self, build_all, tmpdir_factory):
-        filepath = tmpdir_factory.mktemp("custom_math").join("custom-math.tex")
-        build_all.math_documentation.write(filename=filepath)
-        assert Path(filepath).exists()
-
-    @pytest.mark.parametrize(
-        ("filepath", "format"),
-        [(None, "foo"), ("myfile.foo", None), ("myfile.tex", "foo")],
-    )
-    def test_invalid_format(self, build_all, tmpdir_factory, filepath, format):
-        if filepath is not None:
-            filepath = tmpdir_factory.mktemp("custom_math").join(filepath)
-        with pytest.raises(ValueError) as excinfo:  # noqa: PT011
-            build_all.math_documentation.write(filename="foo", format=format)
-        assert check_error_or_warning(
-            excinfo, "Math documentation format must be one of"
-        )
+from .common.util import check_error_or_warning
 
 
 class TestLatexBackendModel:
@@ -465,8 +406,13 @@ def test_create_obj_list(self, dummy_latex_backend_model):
             ),
         ],
     )
-    def test_generate_math_doc(self, dummy_model_data, format, expected):
-        backend_model = latex_backend_model.LatexBackendModel(dummy_model_data)
+    def test_generate_math_doc(
+        self, dummy_model_data, dummy_model_math, format, expected
+    ):
+        backend_model = latex_backend_model.LatexBackendModel(
+            dummy_model_data, dummy_model_math
+        )
+        backend_model._add_all_inputs_as_parameters()
         backend_model.add_global_expression(
             "expr",
             {
@@ -478,8 +424,10 @@ def test_generate_math_doc(self, dummy_model_data, format, expected):
         doc = backend_model.generate_math_doc(format=format)
         assert doc == expected
 
-    def test_generate_math_doc_no_params(self, dummy_model_data):
-        backend_model = latex_backend_model.LatexBackendModel(dummy_model_data)
+    def test_generate_math_doc_no_params(self, dummy_model_data, dummy_model_math):
+        backend_model = latex_backend_model.LatexBackendModel(
+            dummy_model_data, dummy_model_math
+        )
         backend_model.add_global_expression(
             "expr",
             {
@@ -508,8 +456,12 @@ def test_generate_math_doc_no_params(self, dummy_model_data):
                     """
         )
 
-    def test_generate_math_doc_mkdocs_features_tabs(self, dummy_model_data):
-        backend_model = latex_backend_model.LatexBackendModel(dummy_model_data)
+    def test_generate_math_doc_mkdocs_features_tabs(
+        self, dummy_model_data, dummy_model_math
+    ):
+        backend_model = latex_backend_model.LatexBackendModel(
+            dummy_model_data, dummy_model_math
+        )
         backend_model.add_global_expression(
             "expr",
             {
@@ -547,8 +499,13 @@ def test_generate_math_doc_mkdocs_features_tabs(self, dummy_model_data):
                     """
         )
 
-    def test_generate_math_doc_mkdocs_features_admonition(self, dummy_model_data):
-        backend_model = latex_backend_model.LatexBackendModel(dummy_model_data)
+    def test_generate_math_doc_mkdocs_features_admonition(
+        self, dummy_model_data, dummy_model_math
+    ):
+        backend_model = latex_backend_model.LatexBackendModel(
+            dummy_model_data, dummy_model_math
+        )
+        backend_model._add_all_inputs_as_parameters()
         backend_model.add_global_expression(
             "expr",
             {
@@ -600,8 +557,12 @@ def test_generate_math_doc_mkdocs_features_admonition(self, dummy_model_data):
                     """
         )
 
-    def test_generate_math_doc_mkdocs_features_not_in_md(self, dummy_model_data):
-        backend_model = latex_backend_model.LatexBackendModel(dummy_model_data)
+    def test_generate_math_doc_mkdocs_features_not_in_md(
+        self, dummy_model_data, dummy_model_math
+    ):
+        backend_model = latex_backend_model.LatexBackendModel(
+            dummy_model_data, dummy_model_math
+        )
         with pytest.raises(exceptions.ModelError) as excinfo:
             backend_model.generate_math_doc(format="rst", mkdocs_features=True)
 
@@ -718,8 +679,11 @@ def test_get_variable_bounds_string(self, dummy_latex_backend_model):
         }
         assert refs == {"multi_dim_var"}
 
-    def test_param_type(self, dummy_model_data):
-        backend_model = latex_backend_model.LatexBackendModel(dummy_model_data)
+    def test_param_type(self, dummy_model_data, dummy_model_math):
+        backend_model = latex_backend_model.LatexBackendModel(
+            dummy_model_data, dummy_model_math
+        )
+        backend_model._add_all_inputs_as_parameters()
         backend_model.add_global_expression(
             "expr",
             {
diff --git a/tests/test_backend_module.py b/tests/test_backend_module.py
index 92f87b1b..f220a3b9 100644
--- a/tests/test_backend_module.py
+++ b/tests/test_backend_module.py
@@ -7,15 +7,19 @@
 from calliope.exceptions import BackendError
 
 
-@pytest.mark.parametrize("valid_backend", backend.MODEL_BACKENDS)
+@pytest.mark.parametrize("valid_backend", ["pyomo", "gurobi"])
 def test_valid_model_backend(simple_supply, valid_backend):
     """Requesting a valid model backend must result in a backend instance."""
-    backend_obj = backend.get_model_backend(valid_backend, simple_supply._model_data)
+    backend_obj = backend.get_model_backend(
+        valid_backend, simple_supply._model_data, simple_supply.applied_math
+    )
     assert isinstance(backend_obj, BackendModel)
 
 
 @pytest.mark.parametrize("spam", ["not_real", None, True, 1])
-def test_invalid_model_backend(spam):
+def test_invalid_model_backend(spam, simple_supply):
     """Backend requests should catch invalid setups."""
     with pytest.raises(BackendError):
-        backend.get_model_backend(spam, None)
+        backend.get_model_backend(
+            spam, simple_supply._model_data, simple_supply.applied_math
+        )
diff --git a/tests/test_backend_parsing.py b/tests/test_backend_parsing.py
index 5d91f108..8847738c 100644
--- a/tests/test_backend_parsing.py
+++ b/tests/test_backend_parsing.py
@@ -219,14 +219,14 @@ def _equation_slice_obj(name):
 
 
 @pytest.fixture
-def dummy_backend_interface(dummy_model_data):
+def dummy_backend_interface(dummy_model_data, dummy_model_math):
     # ignore the need to define the abstract methods from backend_model.BackendModel
     with patch.multiple(backend_model.BackendModel, __abstractmethods__=set()):
 
         class DummyBackendModel(backend_model.BackendModel):
             def __init__(self):
                 backend_model.BackendModel.__init__(
-                    self, dummy_model_data, instance=None
+                    self, dummy_model_data, dummy_model_math, instance=None
                 )
 
                 self._dataset = dummy_model_data.copy(deep=True)
diff --git a/tests/test_backend_pyomo.py b/tests/test_backend_pyomo.py
index 88f1e0ee..df8d6431 100755
--- a/tests/test_backend_pyomo.py
+++ b/tests/test_backend_pyomo.py
@@ -1,6 +1,4 @@
-import importlib
 import logging
-from copy import deepcopy
 from itertools import product
 
 import numpy as np
@@ -11,9 +9,10 @@
 from pyomo.core.kernel.piecewise_library.transforms import piecewise_sos2
 
 import calliope
+import calliope.backend
 import calliope.exceptions as exceptions
-from calliope.attrdict import AttrDict
-from calliope.backend.pyomo_backend_model import PyomoBackendModel
+import calliope.preprocess
+from calliope.backend import PyomoBackendModel
 
 from .common.util import build_test_model as build_model
 from .common.util import check_error_or_warning, check_variable_exists
@@ -1520,7 +1519,7 @@ def cluster_model(
         override = {
             "config.init.time_subset": ["2005-01-01", "2005-01-04"],
             "config.init.time_cluster": "data_tables/cluster_days.csv",
-            "config.init.add_math": (
+            "config.build.add_math": (
                 ["storage_inter_cluster"] if storage_inter_cluster else []
             ),
             "config.build.cyclic_storage": cyclic,
@@ -1626,59 +1625,34 @@ def simple_supply_updated_cost_flow_cap(
     def temp_path(self, tmpdir_factory):
         return tmpdir_factory.mktemp("custom_math")
 
-    @pytest.mark.parametrize("mode", ["operate", "spores"])
+    @pytest.mark.parametrize("mode", ["operate", "plan"])
     def test_add_run_mode_custom_math(self, caplog, mode):
         caplog.set_level(logging.DEBUG)
-        mode_custom_math = AttrDict.from_yaml(
-            importlib.resources.files("calliope") / "math" / f"{mode}.yaml"
-        )
         m = build_model({}, "simple_supply,two_hours,investment_costs")
+        math = calliope.preprocess.CalliopeMath([mode])
 
-        base_math = deepcopy(m.math)
-        base_math.union(mode_custom_math, allow_override=True)
-
-        backend = PyomoBackendModel(m.inputs, mode=mode)
-        backend._add_run_mode_math()
-
-        assert f"Updating math formulation with {mode} mode math." in caplog.text
+        backend = PyomoBackendModel(m.inputs, math, mode=mode)
 
-        assert m.math != base_math
-        assert backend.inputs.attrs["math"].as_dict() == base_math.as_dict()
+        assert backend.math == math
 
-    def test_add_run_mode_custom_math_before_build(self, caplog, temp_path):
-        """A user can override the run mode math by including it directly in the additional math list"""
+    def test_add_run_mode_custom_math_before_build(self, caplog):
+        """Run mode math is applied before anything else."""
         caplog.set_level(logging.DEBUG)
-        custom_math = AttrDict({"variables": {"flow_cap": {"active": True}}})
-        file_path = temp_path.join("custom-math.yaml")
-        custom_math.to_yaml(file_path)
+        custom_math = {"constraints": {"force_zero_area_use": {"active": True}}}
 
         m = build_model(
-            {"config.init.add_math": ["operate", str(file_path)]},
+            {
+                "config.build.operate_window": "12H",
+                "config.build.operate_horizon": "12H",
+            },
             "simple_supply,two_hours,investment_costs",
         )
-        backend = PyomoBackendModel(m.inputs, mode="operate")
-        backend._add_run_mode_math()
-
-        # We set operate mode explicitly in our additional math so it won't be added again
-        assert "Updating math formulation with operate mode math." not in caplog.text
+        m.build(mode="operate", add_math_dict=custom_math)
 
         # operate mode set it to false, then our math set it back to active
-        assert m.math.variables.flow_cap.active
+        assert m.applied_math.data.constraints.force_zero_area_use.active
         # operate mode set it to false and our math did not override that
-        assert not m.math.variables.storage_cap.active
-
-    def test_run_mode_mismatch(self):
-        m = build_model(
-            {"config.init.add_math": ["operate"]},
-            "simple_supply,two_hours,investment_costs",
-        )
-        backend = PyomoBackendModel(m.inputs)
-        with pytest.warns(exceptions.ModelWarning) as excinfo:
-            backend._add_run_mode_math()
-
-        assert check_error_or_warning(
-            excinfo, "Running in plan mode, but run mode(s) {'operate'}"
-        )
+        assert not m.applied_math.data.variables.storage_cap.active
 
     def test_new_build_get_variable(self, simple_supply):
         """Check a decision variable has the correct data type and has all expected attributes."""
@@ -2255,3 +2229,84 @@ def test_yaml_with_invalid_constraint(self, simple_supply_yaml_invalid):
         )
         # Since we listed only one (invalid) constraint, tracking should not be active
         assert not m.backend.shadow_prices.is_active
+
+
+class TestValidateMathDict:
+    LOGGER = "calliope.backend.backend_model"
+
+    @pytest.fixture
+    def validate_math(self):
+        def _validate_math(math_dict: dict):
+            m = build_model({}, "simple_supply,investment_costs")
+            math = calliope.preprocess.CalliopeMath(["plan", math_dict])
+            backend = calliope.backend.PyomoBackendModel(m._model_data, math)
+            backend._add_all_inputs_as_parameters()
+            backend._validate_math_string_parsing()
+
+        return _validate_math
+
+    def test_base_math(self, caplog, validate_math):
+        with caplog.at_level(logging.INFO, logger=self.LOGGER):
+            validate_math({})
+        assert "Optimisation Model | Validated math strings." in [
+            rec.message for rec in caplog.records
+        ]
+
+    @pytest.mark.parametrize(
+        ("equation", "where"),
+        [
+            ("1 == 1", "True"),
+            (
+                "sum(flow_out * flow_out_eff, over=[nodes, carriers, techs, timesteps]) <= .inf",
+                "base_tech=supply and flow_out_eff>0",
+            ),
+        ],
+    )
+    def test_add_math(self, caplog, validate_math, equation, where):
+        with caplog.at_level(logging.INFO, logger=self.LOGGER):
+            validate_math(
+                {
+                    "constraints": {
+                        "foo": {"equations": [{"expression": equation}], "where": where}
+                    }
+                }
+            )
+        assert "Optimisation Model | Validated math strings." in [
+            rec.message for rec in caplog.records
+        ]
+
+    @pytest.mark.parametrize(
+        "component_dict",
+        [
+            {"equations": [{"expression": "1 = 1"}]},
+            {"equations": [{"expression": "1 = 1"}], "where": "foo[bar]"},
+        ],
+    )
+    @pytest.mark.parametrize("both_fail", [True, False])
+    def test_add_math_fails(self, validate_math, component_dict, both_fail):
+        math_dict = {"constraints": {"foo": component_dict}}
+        errors_to_check = [
+            "math string parsing (marker indicates where parsing stopped, but may not point to the root cause of the issue)",
+            " * constraints:foo:",
+            "equations[0].expression",
+            "where",
+        ]
+        if both_fail:
+            math_dict["constraints"]["bar"] = component_dict
+            errors_to_check.append("* constraints:bar:")
+        else:
+            math_dict["constraints"]["bar"] = {"equations": [{"expression": "1 == 1"}]}
+
+        with pytest.raises(calliope.exceptions.ModelError) as excinfo:
+            validate_math(math_dict)
+        assert check_error_or_warning(excinfo, errors_to_check)
+
+    @pytest.mark.parametrize("eq_string", ["1 = 1", "1 ==\n1[a]"])
+    def test_add_math_fails_marker_correct_position(self, validate_math, eq_string):
+        math_dict = {"constraints": {"foo": {"equations": [{"expression": eq_string}]}}}
+
+        with pytest.raises(calliope.exceptions.ModelError) as excinfo:
+            validate_math(math_dict)
+        errorstrings = str(excinfo.value).split("\n")
+        # marker should be at the "=" sign, i.e., 2 characters from the end
+        assert len(errorstrings[-2]) - 2 == len(errorstrings[-1])
diff --git a/tests/test_core_model.py b/tests/test_core_model.py
index 2d3ad2b1..e16ebfa4 100644
--- a/tests/test_core_model.py
+++ b/tests/test_core_model.py
@@ -1,11 +1,12 @@
 import logging
 from contextlib import contextmanager
 
-import numpy as np
 import pandas as pd
 import pytest
 
 import calliope
+import calliope.backend
+import calliope.preprocess
 
 from .common.util import build_test_model as build_model
 from .common.util import check_error_or_warning
@@ -66,222 +67,6 @@ def test_add_observed_dict_not_dict(self, national_scale_example):
         )
 
 
-class TestAddMath:
-    @pytest.fixture(scope="class")
-    def storage_inter_cluster(self):
-        return build_model(
-            {"config.init.add_math": ["storage_inter_cluster"]},
-            "simple_supply,two_hours,investment_costs",
-        )
-
-    @pytest.fixture(scope="class")
-    def storage_inter_cluster_plus_user_def(self, temp_path, dummy_int: int):
-        new_constraint = calliope.AttrDict(
-            {"variables": {"storage": {"bounds": {"min": dummy_int}}}}
-        )
-        file_path = temp_path.join("custom-math.yaml")
-        new_constraint.to_yaml(file_path)
-        return build_model(
-            {"config.init.add_math": ["storage_inter_cluster", str(file_path)]},
-            "simple_supply,two_hours,investment_costs",
-        )
-
-    @pytest.fixture(scope="class")
-    def temp_path(self, tmpdir_factory):
-        return tmpdir_factory.mktemp("custom_math")
-
-    def test_internal_override(self, storage_inter_cluster):
-        assert "storage_intra_max" in storage_inter_cluster.math["constraints"].keys()
-
-    def test_variable_bound(self, storage_inter_cluster):
-        assert (
-            storage_inter_cluster.math["variables"]["storage"]["bounds"]["min"]
-            == -np.inf
-        )
-
-    @pytest.mark.parametrize(
-        ("override", "expected"),
-        [
-            (["foo"], ["foo"]),
-            (["bar", "foo"], ["bar", "foo"]),
-            (["foo", "storage_inter_cluster"], ["foo"]),
-            (["foo.yaml"], ["foo.yaml"]),
-        ],
-    )
-    def test_allowed_internal_constraint(self, override, expected):
-        with pytest.raises(calliope.exceptions.ModelError) as excinfo:
-            build_model(
-                {"config.init.add_math": override},
-                "simple_supply,two_hours,investment_costs",
-            )
-        assert check_error_or_warning(
-            excinfo,
-            f"Attempted to load additional math that does not exist: {expected}",
-        )
-
-    def test_internal_override_from_yaml(self, temp_path):
-        new_constraint = calliope.AttrDict(
-            {
-                "constraints": {
-                    "constraint_name": {
-                        "foreach": [],
-                        "where": "",
-                        "equations": [{"expression": ""}],
-                    }
-                }
-            }
-        )
-        new_constraint.to_yaml(temp_path.join("custom-math.yaml"))
-        m = build_model(
-            {"config.init.add_math": [str(temp_path.join("custom-math.yaml"))]},
-            "simple_supply,two_hours,investment_costs",
-        )
-        assert "constraint_name" in m.math["constraints"].keys()
-
-    def test_override_existing_internal_constraint(self, temp_path, simple_supply):
-        file_path = temp_path.join("custom-math.yaml")
-        new_constraint = calliope.AttrDict(
-            {
-                "constraints": {
-                    "flow_capacity_per_storage_capacity_min": {"foreach": ["nodes"]}
-                }
-            }
-        )
-        new_constraint.to_yaml(file_path)
-        m = build_model(
-            {"config.init.add_math": [str(file_path)]},
-            "simple_supply,two_hours,investment_costs",
-        )
-        base = simple_supply.math["constraints"][
-            "flow_capacity_per_storage_capacity_min"
-        ]
-        new = m.math["constraints"]["flow_capacity_per_storage_capacity_min"]
-
-        for i in base.keys():
-            if i == "foreach":
-                assert new[i] == ["nodes"]
-            else:
-                assert base[i] == new[i]
-
-    def test_override_order(self, temp_path, simple_supply):
-        to_add = []
-        for path_suffix, foreach in [(1, "nodes"), (2, "techs")]:
-            constr = calliope.AttrDict(
-                {
-                    "constraints.flow_capacity_per_storage_capacity_min.foreach": [
-                        foreach
-                    ]
-                }
-            )
-            filepath = temp_path.join(f"custom-math-{path_suffix}.yaml")
-            constr.to_yaml(filepath)
-            to_add.append(str(filepath))
-
-        m = build_model(
-            {"config.init.add_math": to_add}, "simple_supply,two_hours,investment_costs"
-        )
-
-        base = simple_supply.math["constraints"][
-            "flow_capacity_per_storage_capacity_min"
-        ]
-        new = m.math["constraints"]["flow_capacity_per_storage_capacity_min"]
-
-        for i in base.keys():
-            if i == "foreach":
-                assert new[i] == ["techs"]
-            else:
-                assert base[i] == new[i]
-
-    def test_override_existing_internal_constraint_merge(
-        self, simple_supply, storage_inter_cluster, storage_inter_cluster_plus_user_def
-    ):
-        storage_inter_cluster_math = storage_inter_cluster.math["variables"]["storage"]
-        base_math = simple_supply.math["variables"]["storage"]
-        new_math = storage_inter_cluster_plus_user_def.math["variables"]["storage"]
-        expected = {
-            "title": storage_inter_cluster_math["title"],
-            "description": storage_inter_cluster_math["description"],
-            "default": base_math["default"],
-            "unit": base_math["unit"],
-            "foreach": base_math["foreach"],
-            "where": base_math["where"],
-            "bounds": {
-                "min": new_math["bounds"]["min"],
-                "max": base_math["bounds"]["max"],
-            },
-        }
-
-        assert new_math == expected
-
-
-class TestValidateMathDict:
-    def test_base_math(self, caplog, simple_supply):
-        with caplog.at_level(logging.INFO, logger=LOGGER):
-            simple_supply.validate_math_strings(simple_supply.math)
-        assert "Model: validated math strings" in [
-            rec.message for rec in caplog.records
-        ]
-
-    @pytest.mark.parametrize(
-        ("equation", "where"),
-        [
-            ("1 == 1", "True"),
-            (
-                "flow_out * flow_out_eff + sum(cost, over=costs) <= .inf",
-                "base_tech=supply and flow_out_eff>0",
-            ),
-        ],
-    )
-    def test_add_math(self, caplog, simple_supply, equation, where):
-        with caplog.at_level(logging.INFO, logger=LOGGER):
-            simple_supply.validate_math_strings(
-                {
-                    "constraints": {
-                        "foo": {"equations": [{"expression": equation}], "where": where}
-                    }
-                }
-            )
-        assert "Model: validated math strings" in [
-            rec.message for rec in caplog.records
-        ]
-
-    @pytest.mark.parametrize(
-        "component_dict",
-        [
-            {"equations": [{"expression": "1 = 1"}]},
-            {"equations": [{"expression": "1 = 1"}], "where": "foo[bar]"},
-        ],
-    )
-    @pytest.mark.parametrize("both_fail", [True, False])
-    def test_add_math_fails(self, simple_supply, component_dict, both_fail):
-        math_dict = {"constraints": {"foo": component_dict}}
-        errors_to_check = [
-            "math string parsing (marker indicates where parsing stopped, which might not be the root cause of the issue; sorry...)",
-            " * constraints:foo:",
-            "equations[0].expression",
-            "where",
-        ]
-        if both_fail:
-            math_dict["constraints"]["bar"] = component_dict
-            errors_to_check.append("* constraints:bar:")
-        else:
-            math_dict["constraints"]["bar"] = {"equations": [{"expression": "1 == 1"}]}
-
-        with pytest.raises(calliope.exceptions.ModelError) as excinfo:
-            simple_supply.validate_math_strings(math_dict)
-        assert check_error_or_warning(excinfo, errors_to_check)
-
-    @pytest.mark.parametrize("eq_string", ["1 = 1", "1 ==\n1[a]"])
-    def test_add_math_fails_marker_correct_position(self, simple_supply, eq_string):
-        math_dict = {"constraints": {"foo": {"equations": [{"expression": eq_string}]}}}
-
-        with pytest.raises(calliope.exceptions.ModelError) as excinfo:
-            simple_supply.validate_math_strings(math_dict)
-        errorstrings = str(excinfo.value).split("\n")
-        # marker should be at the "=" sign, i.e., 2 characters from the end
-        assert len(errorstrings[-2]) - 2 == len(errorstrings[-1])
-
-
 class TestOperateMode:
     @contextmanager
     def caplog_session(self, request):
@@ -400,6 +185,34 @@ def test_build_operate_not_allowed_build(self):
             m.build(mode="operate")
 
 
+class TestBuild:
+    @pytest.fixture(scope="class")
+    def init_model(self):
+        return build_model({}, "simple_supply,two_hours,investment_costs")
+
+    def test_ignore_mode_math(self, init_model):
+        init_model.build(ignore_mode_math=True, force=True)
+        assert all(
+            var.obj_type == "parameters"
+            for var in init_model.backend._dataset.data_vars.values()
+        )
+
+    def test_add_math_dict_with_mode_math(self, init_model):
+        init_model.build(
+            add_math_dict={"constraints": {"system_balance": {"active": False}}},
+            force=True,
+        )
+        assert len(init_model.backend.constraints) > 0
+        assert "system_balance" not in init_model.backend.constraints
+
+    def test_add_math_dict_ignore_mode_math(self, init_model):
+        new_var = {
+            "variables": {"foo": {"active": True, "bounds": {"min": -1, "max": 1}}}
+        }
+        init_model.build(add_math_dict=new_var, ignore_mode_math=True, force=True)
+        assert set(init_model.backend.variables) == {"foo"}
+
+
 class TestSolve:
     def test_solve_before_build(self):
         m = build_model({}, "simple_supply,two_hours,investment_costs")
diff --git a/tests/test_core_preprocess.py b/tests/test_core_preprocess.py
index 54592bd5..b0f286f4 100644
--- a/tests/test_core_preprocess.py
+++ b/tests/test_core_preprocess.py
@@ -37,18 +37,18 @@ def test_model_from_dict(self, data_source_dir):
     @pytest.mark.filterwarnings(
         "ignore:(?s).*(links, test_link_a_b_elec) | Deactivated:calliope.exceptions.ModelWarning"
     )
-    def test_valid_scenarios(self):
+    def test_valid_scenarios(self, dummy_int):
         """Test that valid scenario definition from overrides raises no error and results in applied scenario."""
         override = AttrDict.from_yaml_string(
-            """
+            f"""
             scenarios:
                 scenario_1: ['one', 'two']
 
             overrides:
                 one:
-                    techs.test_supply_gas.flow_cap_max: 20
+                    techs.test_supply_gas.flow_cap_max: {dummy_int}
                 two:
-                    techs.test_supply_elec.flow_cap_max: 20
+                    techs.test_supply_elec.flow_cap_max: {dummy_int/2}
 
             nodes:
                 a:
@@ -60,24 +60,29 @@ def test_valid_scenarios(self):
         )
         model = build_model(override_dict=override, scenario="scenario_1")
 
-        assert model._model_def_dict.techs.test_supply_gas.flow_cap_max == 20
-        assert model._model_def_dict.techs.test_supply_elec.flow_cap_max == 20
+        assert (
+            model._model_data.sel(techs="test_supply_gas")["flow_cap_max"] == dummy_int
+        )
+        assert (
+            model._model_data.sel(techs="test_supply_elec")["flow_cap_max"]
+            == dummy_int / 2
+        )
 
-    def test_valid_scenario_of_scenarios(self):
+    def test_valid_scenario_of_scenarios(self, dummy_int):
         """Test that valid scenario definition which groups scenarios and overrides raises
         no error and results in applied scenario.
         """
         override = AttrDict.from_yaml_string(
-            """
+            f"""
             scenarios:
                 scenario_1: ['one', 'two']
                 scenario_2: ['scenario_1', 'new_location']
 
             overrides:
                 one:
-                    techs.test_supply_gas.flow_cap_max: 20
+                    techs.test_supply_gas.flow_cap_max: {dummy_int}
                 two:
-                    techs.test_supply_elec.flow_cap_max: 20
+                    techs.test_supply_elec.flow_cap_max: {dummy_int/2}
                 new_location:
                     nodes.b.techs:
                         test_supply_elec:
@@ -92,8 +97,13 @@ def test_valid_scenario_of_scenarios(self):
         )
         model = build_model(override_dict=override, scenario="scenario_2")
 
-        assert model._model_def_dict.techs.test_supply_gas.flow_cap_max == 20
-        assert model._model_def_dict.techs.test_supply_elec.flow_cap_max == 20
+        assert (
+            model._model_data.sel(techs="test_supply_gas")["flow_cap_max"] == dummy_int
+        )
+        assert (
+            model._model_data.sel(techs="test_supply_elec")["flow_cap_max"]
+            == dummy_int / 2
+        )
 
     def test_invalid_scenarios_dict(self):
         """Test that invalid scenario definition raises appropriate error"""
diff --git a/tests/test_core_util.py b/tests/test_core_util.py
index fb740e4a..8e9175ba 100644
--- a/tests/test_core_util.py
+++ b/tests/test_core_util.py
@@ -184,7 +184,7 @@ def test_invalid_dict(self, to_validate, expected_path):
     @pytest.fixture
     def base_math(self):
         return calliope.AttrDict.from_yaml(
-            Path(calliope.__file__).parent / "math" / "base.yaml"
+            Path(calliope.__file__).parent / "math" / "plan.yaml"
         )
 
     @pytest.mark.parametrize(
@@ -195,7 +195,8 @@ def test_validate_math(self, base_math, dict_path):
             Path(calliope.__file__).parent / "config" / "math_schema.yaml"
         )
         to_validate = base_math.union(
-            calliope.AttrDict.from_yaml(dict_path), allow_override=True
+            calliope.AttrDict.from_yaml(dict_path, allow_override=True),
+            allow_override=True,
         )
         schema.validate_dict(to_validate, math_schema, "")
 
diff --git a/tests/test_io.py b/tests/test_io.py
index 60fb0c5e..b496db6b 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -105,7 +105,7 @@ def test_serialised_list_popped(self, request, serialised_list, model_name):
             ("serialised_nones", ["foo_none", "scenario"]),
             (
                 "serialised_dicts",
-                ["foo_dict", "foo_attrdict", "defaults", "config", "math"],
+                ["foo_dict", "foo_attrdict", "defaults", "config", "applied_math"],
             ),
             ("serialised_sets", ["foo_set", "foo_set_1_item"]),
             ("serialised_single_element_list", ["foo_list_1_item", "foo_set_1_item"]),
@@ -182,7 +182,7 @@ def test_save_csv_not_optimal(self):
             with pytest.warns(exceptions.ModelWarning):
                 model.to_csv(out_path, dropna=False)
 
-    @pytest.mark.parametrize("attr", ["config", "math"])
+    @pytest.mark.parametrize("attr", ["config"])
     def test_dicts_as_model_attrs_and_property(self, model_from_file, attr):
         assert attr in model_from_file._model_data.attrs.keys()
         assert hasattr(model_from_file, attr)
@@ -200,11 +200,9 @@ def test_save_read_solve_save_netcdf(self, model, tmpdir_factory):
         model.to_netcdf(out_path)
         model_from_disk = calliope.read_netcdf(out_path)
 
-        # Ensure _model_def_dict doesn't exist to simulate a re-run via the backend
-        delattr(model_from_disk, "_model_def_dict")
+        # Simulate a re-run via the backend
         model_from_disk.build()
         model_from_disk.solve(force=True)
-        assert not hasattr(model_from_disk, "_model_def_dict")
 
         with tempfile.TemporaryDirectory() as tempdir:
             out_path = os.path.join(tempdir, "model.nc")
diff --git a/tests/test_math.py b/tests/test_math.py
index f27fe107..7887b8a1 100644
--- a/tests/test_math.py
+++ b/tests/test_math.py
@@ -11,6 +11,7 @@
 from .common.util import build_lp, build_test_model
 
 CALLIOPE_DIR: Path = importlib.resources.files("calliope")
+PLAN_MATH: AttrDict = AttrDict.from_yaml(CALLIOPE_DIR / "math" / "plan.yaml")
 
 
 @pytest.fixture(scope="class")
@@ -46,7 +47,7 @@ class TestBaseMath:
 
     @pytest.fixture(scope="class")
     def base_math(self):
-        return AttrDict.from_yaml(CALLIOPE_DIR / "math" / "base.yaml")
+        return AttrDict.from_yaml(CALLIOPE_DIR / "math" / "plan.yaml")
 
     def test_flow_cap(self, compare_lps):
         self.TEST_REGISTER.add("variables.flow_cap")
@@ -79,7 +80,7 @@ def test_storage_max(self, compare_lps):
         self.TEST_REGISTER.add("constraints.storage_max")
         model = build_test_model(scenario="simple_storage,two_hours,investment_costs")
         custom_math = {
-            "constraints": {"storage_max": model.math.constraints.storage_max}
+            "constraints": {"storage_max": PLAN_MATH.constraints.storage_max}
         }
         compare_lps(model, custom_math, "storage_max")
 
@@ -94,7 +95,7 @@ def test_flow_out_max(self, compare_lps):
         )
 
         custom_math = {
-            "constraints": {"flow_out_max": model.math.constraints.flow_out_max}
+            "constraints": {"flow_out_max": PLAN_MATH.constraints.flow_out_max}
         }
         compare_lps(model, custom_math, "flow_out_max")
 
@@ -106,7 +107,7 @@ def test_balance_conversion(self, compare_lps):
         )
         custom_math = {
             "constraints": {
-                "balance_conversion": model.math.constraints.balance_conversion
+                "balance_conversion": PLAN_MATH.constraints.balance_conversion
             }
         }
 
@@ -118,7 +119,7 @@ def test_source_max(self, compare_lps):
             {}, "simple_supply_plus,resample_two_days,investment_costs"
         )
         custom_math = {
-            "constraints": {"my_constraint": model.math.constraints.source_max}
+            "constraints": {"my_constraint": PLAN_MATH.constraints.source_max}
         }
         compare_lps(model, custom_math, "source_max")
 
@@ -129,9 +130,7 @@ def test_balance_transmission(self, compare_lps):
             {"techs.test_link_a_b_elec.one_way": True}, "simple_conversion,two_hours"
         )
         custom_math = {
-            "constraints": {
-                "my_constraint": model.math.constraints.balance_transmission
-            }
+            "constraints": {"my_constraint": PLAN_MATH.constraints.balance_transmission}
         }
         compare_lps(model, custom_math, "balance_transmission")
 
@@ -146,7 +145,7 @@ def test_balance_storage(self, compare_lps):
             "simple_storage,two_hours",
         )
         custom_math = {
-            "constraints": {"my_constraint": model.math.constraints.balance_storage}
+            "constraints": {"my_constraint": PLAN_MATH.constraints.balance_storage}
         }
         compare_lps(model, custom_math, "balance_storage")
 
@@ -261,7 +260,7 @@ def _build_and_compare(
                 overrides = {}
 
             model = build_test_model(
-                {"config.init.add_math": [abs_filepath], **overrides}, scenario
+                {"config.build.add_math": [abs_filepath], **overrides}, scenario
             )
 
             compare_lps(model, custom_math, filename)
diff --git a/tests/test_postprocess_math_documentation.py b/tests/test_postprocess_math_documentation.py
new file mode 100644
index 00000000..fb2558de
--- /dev/null
+++ b/tests/test_postprocess_math_documentation.py
@@ -0,0 +1,59 @@
+from pathlib import Path
+
+import pytest
+
+from calliope.postprocess.math_documentation import MathDocumentation
+
+from .common.util import build_test_model, check_error_or_warning
+
+
+class TestMathDocumentation:
+    @pytest.fixture(scope="class")
+    def no_build(self):
+        model = build_test_model({}, "simple_supply,two_hours,investment_costs")
+        model.build()
+        return model
+
+    @pytest.fixture(scope="class")
+    def build_all(self):
+        model = build_test_model({}, "simple_supply,two_hours,investment_costs")
+        model.build()
+        return MathDocumentation(model, include="all")
+
+    @pytest.fixture(scope="class")
+    def build_valid(self):
+        model = build_test_model({}, "simple_supply,two_hours,investment_costs")
+        model.build()
+        return MathDocumentation(model, include="valid")
+
+    @pytest.mark.parametrize(
+        ("format", "startswith"),
+        [
+            ("tex", "\n\\documentclass{article}"),
+            ("rst", "\nObjective"),
+            ("md", "\n## Objective"),
+        ],
+    )
+    @pytest.mark.parametrize("include", ["build_all", "build_valid"])
+    def test_string_return(self, request, format, startswith, include):
+        math_documentation = request.getfixturevalue(include)
+        string_math = math_documentation.write(format=format)
+        assert string_math.startswith(startswith)
+
+    def test_to_file(self, build_all, tmpdir_factory):
+        filepath = tmpdir_factory.mktemp("custom_math").join("custom-math.tex")
+        build_all.write(filename=filepath)
+        assert Path(filepath).exists()
+
+    @pytest.mark.parametrize(
+        ("filepath", "format"),
+        [(None, "foo"), ("myfile.foo", None), ("myfile.tex", "foo")],
+    )
+    def test_invalid_format(self, build_all, tmpdir_factory, filepath, format):
+        if filepath is not None:
+            filepath = tmpdir_factory.mktemp("custom_math").join(filepath)
+        with pytest.raises(ValueError) as excinfo:  # noqa: PT011
+            build_all.write(filename="foo", format=format)
+        assert check_error_or_warning(
+            excinfo, "Math documentation format must be one of"
+        )
diff --git a/tests/test_preprocess_model_data.py b/tests/test_preprocess_model_data.py
index 391b797e..f15393e2 100644
--- a/tests/test_preprocess_model_data.py
+++ b/tests/test_preprocess_model_data.py
@@ -8,7 +8,7 @@
 
 from calliope import exceptions
 from calliope.attrdict import AttrDict
-from calliope.preprocess import data_tables, load
+from calliope.preprocess import data_tables, scenarios
 from calliope.preprocess.model_data import ModelDataFactory
 
 from .common.util import build_test_model as build_model
@@ -17,11 +17,12 @@
 
 @pytest.fixture
 def model_def():
-    filepath = Path(__file__).parent / "common" / "test_model" / "model.yaml"
-    model_def_dict, model_def_path, _ = load.load_model_definition(
-        filepath.as_posix(), scenario="simple_supply,empty_tech_node"
+    model_def_path = Path(__file__).parent / "common" / "test_model" / "model.yaml"
+    model_dict = AttrDict.from_yaml(model_def_path)
+    model_def_override, _ = scenarios.load_scenario_overrides(
+        model_dict, scenario="simple_supply,empty_tech_node"
     )
-    return model_def_dict, model_def_path
+    return model_def_override, model_def_path
 
 
 @pytest.fixture
diff --git a/tests/test_preprocess_model_math.py b/tests/test_preprocess_model_math.py
new file mode 100644
index 00000000..46af363e
--- /dev/null
+++ b/tests/test_preprocess_model_math.py
@@ -0,0 +1,184 @@
+"""Test the model math handler."""
+
+import logging
+from copy import deepcopy
+from pathlib import Path
+from random import shuffle
+
+import pytest
+
+import calliope
+from calliope.exceptions import ModelError
+from calliope.preprocess import CalliopeMath
+
+
+@pytest.fixture
+def model_math_default():
+    return CalliopeMath([])
+
+
+@pytest.fixture(scope="module")
+def def_path(tmp_path_factory):
+    return tmp_path_factory.mktemp("test_model_math")
+
+
+@pytest.fixture(scope="module")
+def user_math(dummy_int):
+    new_vars = {"variables": {"storage": {"bounds": {"min": dummy_int}}}}
+    new_constr = {
+        "constraints": {
+            "foobar": {"foreach": [], "where": "", "equations": [{"expression": ""}]}
+        }
+    }
+    return calliope.AttrDict(new_vars | new_constr)
+
+
+@pytest.fixture(scope="module")
+def user_math_path(def_path, user_math):
+    file_path = def_path / "custom-math.yaml"
+    user_math.to_yaml(def_path / file_path)
+    return "custom-math.yaml"
+
+
+@pytest.mark.parametrize("invalid_obj", [1, "foo", {"foo": "bar"}, True, CalliopeMath])
+def test_invalid_eq(model_math_default, invalid_obj):
+    """Comparisons should not work with invalid objects."""
+    assert not model_math_default == invalid_obj
+
+
+@pytest.mark.parametrize("modes", [[], ["storage_inter_cluster"]])
+class TestInit:
+    def test_init_order(self, caplog, modes, model_math_default):
+        """Math should be added in order, keeping defaults."""
+        with caplog.at_level(logging.INFO):
+            model_math = CalliopeMath(modes)
+        assert all(
+            f"Math preprocessing | added file '{i}'." in caplog.messages for i in modes
+        )
+        assert model_math_default.history + modes == model_math.history
+
+    def test_init_order_user_math(
+        self, modes, user_math_path, def_path, model_math_default
+    ):
+        """User math order should be respected."""
+        modes = modes + [user_math_path]
+        shuffle(modes)
+        model_math = CalliopeMath(modes, def_path)
+        assert model_math_default.history + modes == model_math.history
+
+    def test_init_user_math_invalid_relative(self, modes, user_math_path):
+        """Init with user math should fail if model definition path is not given for a relative path."""
+        with pytest.raises(ModelError):
+            CalliopeMath(modes + [user_math_path])
+
+    def test_init_user_math_valid_absolute(self, modes, def_path, user_math_path):
+        """Init with user math should succeed if user math is an absolute path."""
+        abs_path = str((def_path / user_math_path).absolute())
+        model_math = CalliopeMath(modes + [abs_path])
+        assert model_math.in_history(abs_path)
+
+    def test_init_dict(self, modes, user_math_path, def_path):
+        """Math dictionary reload should lead to no alterations."""
+        modes = modes + [user_math_path]
+        shuffle(modes)
+        model_math = CalliopeMath(modes, def_path)
+        saved = dict(model_math)
+        reloaded = CalliopeMath.from_dict(saved)
+        assert model_math == reloaded
+
+
+class TestMathLoading:
+    @pytest.fixture(scope="class")
+    def pre_defined_mode(self):
+        return "storage_inter_cluster"
+
+    @pytest.fixture
+    def model_math_w_mode(self, model_math_default, pre_defined_mode):
+        model_math_default._add_pre_defined_file(pre_defined_mode)
+        return model_math_default
+
+    @pytest.fixture
+    def model_math_w_mode_user(self, model_math_w_mode, user_math_path, def_path):
+        model_math_w_mode._add_user_defined_file(user_math_path, def_path)
+        return model_math_w_mode
+
+    @pytest.fixture(scope="class")
+    def predefined_mode_data(self, pre_defined_mode):
+        path = Path(calliope.__file__).parent / "math" / f"{pre_defined_mode}.yaml"
+        math = calliope.AttrDict.from_yaml(path)
+        return math
+
+    def test_predefined_add(self, model_math_w_mode, predefined_mode_data):
+        """Added mode should be in data."""
+        flat = predefined_mode_data.as_dict_flat()
+        assert all(model_math_w_mode.data.get_key(i) == flat[i] for i in flat.keys())
+
+    def test_predefined_add_history(self, pre_defined_mode, model_math_w_mode):
+        """Added modes should be recorded."""
+        assert model_math_w_mode.in_history(pre_defined_mode)
+
+    def test_predefined_add_duplicate(self, pre_defined_mode, model_math_w_mode):
+        """Adding the same mode twice is invalid."""
+        with pytest.raises(ModelError):
+            model_math_w_mode._add_pre_defined_file(pre_defined_mode)
+
+    @pytest.mark.parametrize("invalid_mode", ["foobar", "foobar.yaml", "operate.yaml"])
+    def test_predefined_add_fail(self, invalid_mode, model_math_w_mode):
+        """Requesting inexistent modes or modes with suffixes should fail."""
+        with pytest.raises(ModelError):
+            model_math_w_mode._add_pre_defined_file(invalid_mode)
+
+    def test_user_math_add(
+        self, model_math_w_mode_user, predefined_mode_data, user_math
+    ):
+        """Added user math should be in data."""
+        expected_math = deepcopy(predefined_mode_data)
+        expected_math.union(user_math, allow_override=True)
+        flat = expected_math.as_dict_flat()
+        assert all(
+            model_math_w_mode_user.data.get_key(i) == flat[i] for i in flat.keys()
+        )
+
+    def test_user_math_add_history(self, model_math_w_mode_user, user_math_path):
+        """Added user math should be recorded."""
+        assert model_math_w_mode_user.in_history(user_math_path)
+
+    def test_repr(self, model_math_w_mode):
+        expected_repr_content = """Calliope math definition dictionary with:
+    4 decision variable(s)
+    0 global expression(s)
+    9 constraint(s)
+    0 piecewise constraint(s)
+    0 objective(s)
+        """
+        assert expected_repr_content == str(model_math_w_mode)
+
+    def test_add_dict(self, model_math_w_mode, model_math_w_mode_user, user_math):
+        model_math_w_mode.add(user_math)
+        assert model_math_w_mode_user == model_math_w_mode
+
+    def test_user_math_add_duplicate(
+        self, model_math_w_mode_user, user_math_path, def_path
+    ):
+        """Adding the same user math file twice should fail."""
+        with pytest.raises(ModelError):
+            model_math_w_mode_user._add_user_defined_file(user_math_path, def_path)
+
+    @pytest.mark.parametrize("invalid_mode", ["foobar", "foobar.yaml", "operate.yaml"])
+    def test_user_math_add_fail(self, invalid_mode, model_math_w_mode_user, def_path):
+        """Requesting inexistent user modes should fail."""
+        with pytest.raises(ModelError):
+            model_math_w_mode_user._add_user_defined_file(invalid_mode, def_path)
+
+
+class TestValidate:
+    def test_validate_math_fail(self):
+        """Invalid math keys must trigger a failure."""
+        model_math = CalliopeMath([{"foo": "bar"}])
+        with pytest.raises(ModelError):
+            model_math.validate()
+
+    def test_math_default(self, caplog, model_math_default):
+        with caplog.at_level(logging.INFO):
+            model_math_default.validate()
+        assert "Math preprocessing | validated math against schema." in caplog.messages
diff --git a/tests/test_preprocess_time.py b/tests/test_preprocess_time.py
index 942d086f..3272872d 100644
--- a/tests/test_preprocess_time.py
+++ b/tests/test_preprocess_time.py
@@ -111,7 +111,7 @@ def test_cluster_datesteps(self, clustered_model):
     @pytest.mark.parametrize(
         "var",
         [
-            "lookup_cluster_first_timestep",
+            "cluster_first_timestep",
             "lookup_cluster_last_timestep",
             "lookup_datestep_cluster",
             "lookup_datestep_last_cluster_timestep",