Repository: calliope-project/calliope (language: python)

Issue:
Add public API to update a built optimisation model
## Problem description

Since we have YAML math we also have a new set of methods to build optimisation models using Pyomo (and hopefully Gurobi python API). Certain methods need to be made available on the model object so a user can do the following:

- [x] Check the content of a specific variable/constraint/expression
- [x] Update a parameter's values (and update the corresponding expressions/constraints/objectives that rely on that parameter)
- [x] Fix a variable, so it acts as a parameter
- [ ] (De)activate a constraint so it is(n't) applied in the optimisation
- [x] Update decision variable bounds

This replaces #293 for the 0.7.x development (#293 will stay open for any future v0.6.x dev).

## Calliope version

v0.7.0-dev



Discussion / hints:

Fix patch:
diff --git a/doc/api/api.rst b/doc/api/api.rst
index 91aed705b..ec2dda34a 100644
--- a/doc/api/api.rst
+++ b/doc/api/api.rst
@@ -33,6 +33,15 @@ Time series
 .. automodule:: calliope.time.funcs
     :members: resample
 
+.. _api_backend_interface:
+
+Methods to interface with the optimisation problem (a.k.a., `backend`)
+======================================================================
+
+.. autoclass:: calliope.backend.backend_model.BackendModel
+    :members:
+    :inherited-members:
+
 .. _api_utility_classes:
 
 Utility classes: AttrDict, Exceptions, Logging
diff --git a/doc/user/troubleshooting.rst b/doc/user/troubleshooting.rst
index 8ffe00226..6691f7a9e 100644
--- a/doc/user/troubleshooting.rst
+++ b/doc/user/troubleshooting.rst
@@ -20,7 +20,8 @@ General strategies
 
   .. code-block:: python
 
-    model.to_lp('my_saved_model.lp')
+    model.build()
+    model.backend.to_lp('my_saved_model.lp')
 
 Improving solution times
 ------------------------
diff --git a/src/calliope/backend/backend_model.py b/src/calliope/backend/backend_model.py
new file mode 100644
index 000000000..5fdd5790f
--- /dev/null
+++ b/src/calliope/backend/backend_model.py
@@ -0,0 +1,767 @@
+# Copyright (C) since 2013 Calliope contributors listed in AUTHORS.
+# Licensed under the Apache 2.0 License (see LICENSE file).
+
+from __future__ import annotations
+
+import logging
+import typing
+from abc import ABC, abstractmethod
+from pathlib import Path
+from typing import (
+    TYPE_CHECKING,
+    Any,
+    Callable,
+    Generic,
+    Literal,
+    Optional,
+    SupportsFloat,
+    TypeVar,
+    Union,
+)
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+
+from calliope.backend import parsing
+
+if TYPE_CHECKING:
+    from calliope.backend.parsing import T as Tp
+
+from calliope.exceptions import BackendError
+
+T = TypeVar("T")
+_COMPONENTS_T = Literal[
+    "variables", "constraints", "objectives", "parameters", "global_expressions"
+]
+
+logger = logging.getLogger(__name__)
+
+
+class BackendModelGenerator(ABC):
+    _VALID_COMPONENTS: tuple[_COMPONENTS_T, ...] = typing.get_args(_COMPONENTS_T)
+    _COMPONENT_ATTR_METADATA = ["description", "unit"]
+
+    def __init__(self, inputs: xr.Dataset):
+        """Abstract base class to build a representation of the optimisation problem.
+
+        Args:
+            inputs (xr.Dataset): Calliope model data.
+        """
+
+        self._dataset = xr.Dataset()
+        self.inputs = inputs
+        self.valid_math_element_names: set = set()
+
+    @abstractmethod
+    def add_parameter(
+        self,
+        parameter_name: str,
+        parameter_values: xr.DataArray,
+        default: Any = np.nan,
+        use_inf_as_na: bool = False,
+    ) -> None:
+        """
+        Add input parameter to backend model in-place.
+        If the backend interface allows for mutable parameter objects, they will be
+        generated, otherwise a copy of the model input dataset will be used.
+        In either case, NaN values are filled with the given parameter default value.
+
+        Args:
+            parameter_name (str): Name of parameter.
+            parameter_values (xr.DataArray): Array of parameter values.
+            default (Any, optional):
+                Default value to fill NaN entries in parameter values array.
+                Defaults to np.nan.
+            use_inf_as_na (bool, optional):
+                If True, will consider np.inf parameter value entries as np.nan and
+                consequently try to fill those entries with the parameter default value.
+                Defaults to False.
+        """
+
+    @abstractmethod
+    def add_constraint(
+        self,
+        name: str,
+        constraint_dict: parsing.UnparsedConstraintDict,
+    ) -> None:
+        """
+        Add constraint equation to backend model in-place.
+        Resulting backend dataset entries will be constraint objects.
+
+        Args:
+            model_data (xr.Dataset):
+                Calliope model data with which to create an array mask - only those
+                dataset entries in the mask will be generated.
+            name (str):
+                Name of the constraint
+            constraint_dict (parsing.UnparsedConstraintDict):
+                Constraint configuration dictionary, ready to be parsed and then evaluated.
+        """
+
+    @abstractmethod
+    def add_global_expression(
+        self,
+        name: str,
+        expression_dict: parsing.UnparsedExpressionDict,
+    ) -> None:
+        """
+        Add global expression (arithmetic combination of parameters and/or decision variables)
+        to backend model in-place.
+        Resulting backend dataset entries will be linear expression objects.
+
+        Args:
+            model_data (xr.Dataset):
+                Calliope model data with which to create an array mask - only those
+                dataset entries in the mask will be generated.
+            name (str):
+                Name of the global expression
+            expression_dict (parsing.UnparsedExpressionDict):
+                Global expression configuration dictionary, ready to be parsed and then evaluated.
+        """
+
+    @abstractmethod
+    def add_variable(
+        self,
+        name: str,
+        variable_dict: parsing.UnparsedVariableDict,
+    ) -> None:
+        """
+        Add decision variable to backend model in-place.
+        Resulting backend dataset entries will be decision variable objects.
+
+        Args:
+            model_data (xr.Dataset):
+                Calliope model data with which to create an array mask - only those
+                dataset entries in the mask will be generated.
+            name (str):
+                Name of the variable.
+            variable_dict (parsing.UnparsedVariableDict):
+                Variable configuration dictionary, ready to be parsed and then evaluated.
+        """
+
+    @abstractmethod
+    def add_objective(
+        self,
+        name: str,
+        objective_dict: parsing.UnparsedObjectiveDict,
+    ) -> None:
+        """
+        Add objective arithmetic to backend model in-place.
+        Resulting backend dataset entry will be a single, unindexed objective object.
+
+        Args:
+            model_data (xr.Dataset):
+                Calliope model data with which to create a constraint mask - only those
+                dataset entries in the mask will be generated.
+            name (str):
+                Name of the objective.
+            objective_dict (parsing.UnparsedObjectiveDict):
+                Objective configuration dictionary, ready to be parsed and then evaluated.
+        """
+
+    def _add_component(
+        self,
+        name: str,
+        component_dict: Optional[Tp],
+        component_setter: Callable,
+        component_type: Literal[
+            "variables", "global_expressions", "constraints", "objectives"
+        ],
+        break_early: bool = True,
+    ) -> Optional[parsing.ParsedBackendComponent]:
+        """Generalised function to add a optimisation problem component array to the model.
+
+        Args:
+            name:
+                Name of the component.
+                If not providing the `component_dict` directly, this name must be available in the input math provided on initialising the class.
+            component_dict (Union[parsing.UnparsedConstraintDict, parsing.UnparsedExpressionDict, parsing.UnparsedObjectiveDict, parsing.UnparsedVariableDict]):
+                Unparsed YAML dictionary configuration.
+            component_setter (Callable):
+                Function to combine evaluated xarray DataArrays into backend component objects.
+                Will receive outputs of `evaluate_where` and (optionally, if `equations` is a component dictionary key) `evaluate_expression` as inputs.
+            component_type (Literal["variables", "global_expressions", "constraints", "objectives"])
+
+        Raises:
+            BackendError:
+                The sub-equations of the parsed component cannot generate component
+                objects on duplicate index entries.
+        """
+        references: set[str] = set()
+
+        if component_dict is None:
+            component_dict = self.inputs.math[component_type][name]
+
+        if component_dict.get("active", False):
+            logger.debug(
+                f"({component_type.removesuffix('s')}, {name}): "
+                "Component deactivated and therefore not built."
+            )
+            return None
+
+        self._raise_error_on_preexistence(name, component_type)
+        parsed_component = parsing.ParsedBackendComponent(
+            component_type, name, component_dict
+        )
+
+        top_level_where = parsed_component.generate_top_level_where_array(
+            self.inputs, align_to_foreach_sets=False, break_early=break_early
+        )
+        if break_early and not top_level_where.any():
+            return parsed_component
+
+        self._create_obj_list(name, component_type)
+
+        equations = parsed_component.parse_equations(self.valid_math_element_names)
+        if not equations:
+            component_da = component_setter(
+                parsed_component.drop_dims_not_in_foreach(top_level_where)
+            )
+        else:
+            component_da = (
+                xr.DataArray()
+                .where(parsed_component.drop_dims_not_in_foreach(top_level_where))
+                .astype(np.dtype("O"))
+            )
+        for element in equations:
+            where = element.evaluate_where(self.inputs, initial_where=top_level_where)
+            if break_early and not where.any():
+                continue
+
+            where = parsed_component.drop_dims_not_in_foreach(where)
+
+            if component_da.where(where).notnull().any():
+                if component_da.shape:
+                    overlap = component_da.where(where).to_series().dropna().index
+                    substring = (
+                        f"trying to set two equations for the same index:\n{overlap}"
+                    )
+                else:
+                    substring = "trying to set two equations for the same component."
+
+                self.delete_component(name, component_type)
+                raise BackendError(
+                    f"({component_type.removesuffix('s')}, {element.name}): {substring}"
+                )
+
+            to_fill = component_setter(element, where, references)
+            component_da = component_da.fillna(to_fill)
+
+        if break_early and component_da.isnull().all():
+            self.delete_component(name, component_type)
+            return parsed_component
+
+        self._add_to_dataset(
+            name, component_da, component_type, component_dict, references
+        )
+
+        return parsed_component
+
+    @abstractmethod
+    def delete_component(self, key: str, component_type: _COMPONENTS_T) -> None:
+        """Delete a list object from the backend model object.
+
+        Args:
+            key (str): Name of object.
+            component_type (str): Object type.
+        """
+
+    @abstractmethod
+    def _create_obj_list(self, key: str, component_type: _COMPONENTS_T) -> None:
+        """Attach an empty list object to the backend model object.
+        This may be a backend-specific subclass of a standard list object.
+
+        Args:
+            key (str): Name of object.
+            component_type (str): Object type.
+
+        Raises:
+            BackendError: Cannot overwrite object of same name and type.
+        """
+
+    def _add_all_inputs_as_parameters(self) -> None:
+        """
+        Add all parameters to backend dataset in-place, including those in the run configuration.
+        If model data does not include a parameter, their default values will be added here
+        as unindexed backend dataset parameters.
+
+        TODO: Move the decision on which run config params to generate as backend params
+              earlier in the process.
+        Parameters in "objective_options" and the bigM parameter will be added from
+        run configuration.
+
+        Args:
+            model_data (xr.Dataset): Input model data.
+            defaults (dict): Parameter defaults.
+            run_config (UpdateObserverDict): Run configuration dictionary.
+        """
+
+        for param_name, param_data in self.inputs.filter_by_attrs(
+            is_result=0
+        ).data_vars.items():
+            default_val = self.inputs.attrs["defaults"].get(param_name, np.nan)
+            self.add_parameter(param_name, param_data, default_val)
+        for param_name, default_val in self.inputs.attrs["defaults"].items():
+            if param_name in self.parameters.keys():
+                continue
+            self.add_parameter(
+                param_name, xr.DataArray(default_val), use_inf_as_na=False
+            )
+
+        for option_name, option_val in self.inputs.run_config.get(
+            "objective_options", {}
+        ).items():
+            if option_name == "cost_class":
+                objective_cost_class = {
+                    k: v for k, v in option_val.items() if k in self.inputs.costs
+                }
+                self.add_parameter(
+                    "objective_cost_class",
+                    xr.DataArray.from_series(
+                        pd.Series(objective_cost_class).rename_axis(index="costs")
+                    ),
+                )
+            else:
+                self.add_parameter("objective_" + option_name, xr.DataArray(option_val))
+        self.add_parameter(
+            "bigM", xr.DataArray(self.inputs.run_config.get("bigM", 1e10))
+        )
+
+    @staticmethod
+    def _clean_arrays(*args) -> None:
+        """
+        Preemptively delete objects with large memory footprints that might otherwise
+        stick around longer than necessary.
+        """
+        del args
+
+    def _add_to_dataset(
+        self,
+        name: str,
+        da: xr.DataArray,
+        obj_type: _COMPONENTS_T,
+        unparsed_dict: Union[parsing.UNPARSED_DICTS, dict],
+        references: Optional[set] = None,
+    ):
+        """
+        Add array of backend objects to backend dataset in-place.
+
+        Args:
+            name (str): Name of entry in dataset.
+            da (xr.DataArray): Data to add.
+            obj_type (str): Type of backend objects in the array.
+            unparsed_dict (DT):
+                Dictionary describing the object being added, from which descriptor attributes will be extracted and added to the array attributes.
+            references (set):
+                All other backend objects which are references in this backend object's linear expression(s).
+                E.g. the constraint "carrier_prod / energy_eff <= energy_cap" references the variables ["carrier_prod", "energy_cap"]
+                and the parameter ["energy_eff"].
+                All referenced objects will have their "references" attribute updated with this object's name.
+                Defaults to None.
+        """
+
+        add_attrs = {
+            attr: unparsed_dict.get(attr)
+            for attr in self._COMPONENT_ATTR_METADATA
+            if attr in unparsed_dict.keys()
+        }
+
+        da.attrs.update(
+            {
+                "obj_type": obj_type,
+                "references": set(),
+                "coords_in_name": False,
+                **add_attrs,  # type: ignore
+            }
+        )
+        self._dataset[name] = da
+
+        if references is not None:
+            for reference in references:
+                self._dataset[reference].attrs["references"].add(name)
+
+    def _apply_func(
+        self, func: Callable, *args, output_core_dims: tuple = ((),), **kwargs
+    ) -> xr.DataArray:
+        """
+        Apply a function to every element of an arbitrary number of xarray DataArrays.
+
+        Args:
+            func (Callable):
+                Un-vectorized function to call.
+                Number of accepted args should equal len(args).
+                Number of accepted kwargs should equal len(kwargs).
+            args (xr.DataArray):
+                xarray DataArrays which will be broadcast together and then iterated over
+                to apply the function.
+            output_core_dims (tuple):
+                Additional dimensions which are expected to be passed back from `xr.apply_ufunc` after applying `func`.
+                This is directly passed to `xr.apply_ufunc`; see their documentation for more details.
+                Defaults to ((), )
+            kwargs (dict[str, Any]):
+                Additional keyword arguments to pass to `func`.
+
+        Returns:
+            xr.DataArray: Array with func applied to all elements.
+        """
+        return xr.apply_ufunc(
+            func,
+            *args,
+            kwargs=kwargs,
+            vectorize=True,
+            keep_attrs=True,
+            output_dtypes=[np.dtype("O")],
+            output_core_dims=output_core_dims,
+        )
+
+    def _raise_error_on_preexistence(self, key: str, obj_type: _COMPONENTS_T):
+        """
+        We do not allow any overlap of backend object names since they all have to
+        co-exist in the backend dataset.
+        I.e., users cannot overwrite any backend component with another
+        (of the same type or otherwise).
+
+        Args:
+            key (str): Backend object name
+            obj_type (Literal["variables", "constraints", "objectives", "parameters", "expressions"]): Object type.
+
+        Raises:
+            BackendError:
+                Raised if `key` already exists in the backend model
+                (either with the same or different type as `obj_type`).
+        """
+        if key in self._dataset.keys():
+            if key in getattr(self, obj_type):
+                raise BackendError(
+                    f"Trying to add already existing `{key}` to backend model {obj_type}."
+                )
+            else:
+                other_obj_type = self._dataset[key].attrs["obj_type"].removesuffix("s")
+                raise BackendError(
+                    f"Trying to add already existing *{other_obj_type}* `{key}` "
+                    f"as a backend model *{obj_type.removesuffix('s')}*."
+                )
+
+    @property
+    def constraints(self):
+        "Slice of backend dataset to show only built constraints"
+        return self._dataset.filter_by_attrs(obj_type="constraints")
+
+    @property
+    def variables(self):
+        "Slice of backend dataset to show only built variables"
+        return self._dataset.filter_by_attrs(obj_type="variables")
+
+    @property
+    def parameters(self):
+        "Slice of backend dataset to show only built parameters"
+        return self._dataset.filter_by_attrs(obj_type="parameters")
+
+    @property
+    def global_expressions(self):
+        "Slice of backend dataset to show only built global expressions"
+        return self._dataset.filter_by_attrs(obj_type="global_expressions")
+
+    @property
+    def objectives(self):
+        "Slice of backend dataset to show only built objectives"
+        return self._dataset.filter_by_attrs(obj_type="objectives")
+
+
+class BackendModel(BackendModelGenerator, Generic[T]):
+    def __init__(self, inputs: xr.Dataset, instance: T) -> None:
+        """Abstract base class to build bakcend models that interface with solvers.
+
+        Args:
+            inputs (xr.Dataset): Calliope model data.
+            instance (T): Interface model instance.
+        """
+        super().__init__(inputs)
+        self._instance = instance
+
+    @abstractmethod
+    def get_parameter(
+        self, parameter_name: str, as_backend_objs: bool = True
+    ) -> xr.DataArray:
+        """
+        Extract parameter from backend dataset.
+
+        Args:
+            parameter_name (str): Name of parameter.
+            TODO: hide this and create a method to edit parameter values
+                  (to handle interfaces with non-mutable params)
+            as_backend_objs (bool, optional):
+                If True, will keep the array entries as backend interface objects,
+                which can be updated to update the underlying model.
+                Otherwise, parameter values are given directly, with default values in place of NaNs.
+                Defaults to True.
+
+        Returns:
+            xr.DataArray: parameter array.
+        """
+
+    @abstractmethod
+    def get_constraint(
+        self,
+        constraint_name: str,
+        as_backend_objs: bool = True,
+        eval_body: bool = False,
+    ) -> Union[xr.DataArray, xr.Dataset]:
+        """
+        Get constraint data as either a table of details or as an array of backend interface objects.
+        Can be used to inspect and debug built constraints.
+
+        Args:
+            constraint_name (str): Name of constraint, as given in YAML constraint key.
+            TODO: hide this and create a method to edit constraints that handles differences in interface APIs.
+            as_backend_objs (bool, optional):
+                If True, will keep the array entries as backend interface objects,
+                which can be updated to change the underlying model.
+                Otherwise, constraint body, and lower and upper bounds are given in a table.
+                Defaults to True.
+            eval_body (bool, optional):
+                If True and as_backend_objs is False, will attempt to evaluate the constraint body.
+                If the model has been optimised, this attempt will produce a numeric value to see where the constraint sits between the lower or upper bound.
+                If the model has not yet been optimised, this attempt will fall back on the same as
+                if `eval_body` was set to False, i.e. a string representation of the linear expression in the constraint body.
+                Defaults to False.
+
+        Returns:
+            Union[xr.DataArray, xr.Dataset]:
+                If as_backend_objs is True, will return an xr.DataArray.
+                Otherwise, a xr.Dataset will be given, indexed over the same dimensions as the xr.DataArray, with variables for the constraint body, and upper (`ub`) and lower (`lb`) bounds.
+        """
+
+    @abstractmethod
+    def get_variable(
+        self, variable_name: str, as_backend_objs: bool = True
+    ) -> xr.DataArray:
+        """Extract decision variable array from backend dataset
+
+        Args:
+            variable_name (str): Name of variable.
+            TODO: hide this and create a method to edit variables that handles differences in interface APIs.
+            as_backend_objs (bool, optional):
+                If True, will keep the array entries as backend interface objects,
+                which can be updated to update the underlying model.
+                Otherwise, variable values are given directly.
+                If the model has not been successfully optimised, variable values will all be None.
+                Defaults to True.
+
+        Returns:
+            xr.DataArray: Decision variable array.
+        """
+
+    @abstractmethod
+    def get_variable_bounds(self, name: str) -> xr.Dataset:
+        """Extract decision variable upper and lower bound array from backend dataset
+
+        Args:
+            variable_name (str): Name of variable.
+
+        Returns:
+            xr.Dataset: Contains the arrays for upper ("ub", a.k.a. "max") and lower ("lb", a.k.a. "min") variable bounds.
+        """
+
+    @abstractmethod
+    def get_global_expression(
+        self, expression_name: str, as_backend_objs: bool = True, eval_body: bool = True
+    ) -> xr.DataArray:
+        """Extract global expression array from backend dataset
+
+        Args:
+            global_expression_name (str): Name of global expression
+            TODO: hide this and create a method to edit expressions that handles differences in interface APIs.
+            as_backend_objs (bool, optional):
+                If True, will keep the array entries as backend interface objects,
+                which can be updated to update the underlying model.
+                Otherwise, global expression values are given directly.
+                If the model has not been successfully optimised, expression values will all be provided as strings.
+                Defaults to True.
+            eval_body (bool, optional):
+                If True and as_backend_objs is False, will attempt to evaluate the expression.
+                If the model has been optimised, this attempt will produce a numeric value.
+                If the model has not yet been optimised, this attempt will fall back on the same as
+                if `eval_body` was set to False, i.e. a string representation of the linear expression.
+                Defaults to True.
+
+        Returns:
+            xr.DataArray: global expression array.
+        """
+
+    @abstractmethod
+    def update_parameter(
+        self, name: str, new_values: Union[xr.DataArray, SupportsFloat]
+    ) -> None:
+        """Update parameter elements using an array of new values.
+        If the parameter has not been previously defined, it will be added to the optimisation problem based on the new values given (with NaNs reverting to default values).
+        If the new values have fewer dimensions than are on the parameter array, the new values will be broadcast across the missing dimensions before applying the update.
+
+        Args:
+            name (str): Parameter to update
+            new_values (Union[xr.DataArray, SupportsFloat]): New values to apply. Any empty (NaN) elements in the array will be skipped.
+        """
+
+    @abstractmethod
+    def update_variable_bounds(
+        self,
+        name: str,
+        *,
+        min: Optional[Union[xr.DataArray, SupportsFloat]] = None,
+        max: Optional[Union[xr.DataArray, SupportsFloat]] = None,
+    ) -> None:
+        """
+        Update the bounds on a decision variable.
+        If the variable bounds are defined by parameters in the math formulation,
+        the parameters themselves will be updated.
+
+        Args:
+            name (str): Variable to update.
+            min (Union[xr.DataArray, SupportsFloat], optional):
+                If provided, the Non-NaN values in the array will be used to defined new lower bounds in the decision variable.
+                Defaults to None.
+            max (Union[xr.DataArray, SupportsFloat], optional):
+                If provided, the Non-NaN values in the array will be used to defined new upper bounds in the decision variable.
+                Defaults to None.
+        """
+
+    @abstractmethod
+    def fix_variable(self, name: str, where: Optional[xr.DataArray] = None) -> None:
+        """
+        Fix the variable value to the value quantified on the most recent call to `solve`.
+        Fixed variables will be treated as parameters in the optimisation.
+
+        Args:
+            name (str): Variable to update.
+            where (xr.DataArray, optional):
+                If provided, only a subset of the coordinates in the variable will be fixed.
+                Must be a boolean array or a float equivalent, where NaN is used instead of False.
+                Defaults to None
+        """
+
+    @abstractmethod
+    def unfix_variable(self, name: str, where: Optional[xr.DataArray] = None) -> None:
+        """
+        Unfix the variable so that it is treated as a decision variable in the next call to `solve`.
+
+        Args:
+            name (str): Variable to update
+            where (xr.DataArray, optional):
+                If provided, only a subset of the coordinates in the variable will be unfixed.
+                Must be a boolean array or a float equivalent, where NaN is used instead of False.
+                Defaults to None
+        """
+
+    @abstractmethod
+    def solve(
+        self,
+        solver: str,
+        solver_io: Optional[str] = None,
+        solver_options: Optional[dict] = None,
+        save_logs: Optional[str] = None,
+        warmstart: bool = False,
+        **solve_kwargs,
+    ):
+        """
+        Optimise built model. If solution is optimal, interface objects
+        (decision variables, global expressions, constraints, objective) can be successfully
+        evaluated for their values at optimality.
+
+        Args:
+            solver (str): Name of solver to optimise with.
+            solver_io (Optional[str], optional):
+                If chosen solver has a python interface, set to "python" for potential
+                performance gains, otherwise should be left as None. Defaults to None.
+            solver_options (Optional[dict], optional):
+                Solver options/parameters to pass directly to solver.
+                See solver documentation for available parameters that can be influenced.
+                Defaults to None.
+            save_logs (Optional[str], optional):
+                If given, solver logs and built LP file will be saved to this filepath.
+                Defaults to None.
+            warmstart (bool, optional):
+                If True, and the chosen solver is capable of implementing it, an existing
+                optimal solution will be used to warmstart the next solve run.
+                Defaults to False.
+        """
+
+    def load_results(self) -> xr.Dataset:
+        """
+        Evaluate backend decision variables, global expressions, and parameters (if not in inputs)
+        after a successful model run.
+
+        Returns:
+            xr.Dataset: Dataset of optimal solution results (all numeric data).
+        """
+
+        def _drop_attrs(da):
+            da.attrs = {
+                k: v for k, v in da.attrs.items() if k in self._COMPONENT_ATTR_METADATA
+            }
+            return da
+
+        all_variables = {
+            name_: _drop_attrs(self.get_variable(name_, as_backend_objs=False))
+            for name_, var in self.variables.items()
+            if var.notnull().any()
+        }
+        all_global_expressions = {
+            name_: _drop_attrs(
+                self.get_global_expression(name_, as_backend_objs=False, eval_body=True)
+            )
+            for name_, expr in self.global_expressions.items()
+            if expr.notnull().any()
+        }
+
+        results = xr.Dataset({**all_variables, **all_global_expressions})
+
+        return results
+
+    @abstractmethod
+    def verbose_strings(self) -> None:
+        """
+        Update optimisation model object string representations to include the index coordinates of the object.
+
+        E.g., `variables(carrier_prod)[0]` will become `variables(carrier_prod)[power, region1, ccgt, 2005-01-01 00:00]`
+
+        This takes approximately 10% of the peak memory required to initially build the optimisation problem, so should only be invoked if inspecting the model in detail (e.g., debugging)
+
+        Only string representations of model parameters and variables will be updated since global expressions automatically show the string representation of their contents.
+        """
+
+    @abstractmethod
+    def to_lp(self, path: Union[str, Path]) -> None:
+        """Write the optimisation problem to file in the linear programming LP format.
+        The LP file can be used for debugging and to submit to solvers directly.
+
+        Args:
+            path (Union[str, Path]): Path to which the LP file will be written.
+        """
+
+    def _find_all_references(self, initial_references: set) -> set:
+        """
+        Find all nested references to optimisation problem components from an initial set of references.
+
+        Args:
+            initial_references (set): names of optimisation problem components.
+
+        Returns:
+            set: `initial_references` + any names of optimisation problem components referenced from those initial references.
+        """
+
+        references = initial_references.copy()
+        for reference in initial_references:
+            new_refs = self._dataset[reference].attrs.get("references", {})
+            references.update(self._find_all_references(new_refs))
+        return references
+
+    def _rebuild_reference(self, reference: str) -> None:
+        """Delete and rebuild an optimisation problem component.
+
+        Args:
+            references (str): name of optimisation problem component.
+        """
+        obj_type = self._dataset[reference].attrs["obj_type"]
+        self.delete_component(reference, obj_type)
+        getattr(self, "add_" + obj_type.removesuffix("s"))(name=reference)
diff --git a/src/calliope/backend/backends.py b/src/calliope/backend/backends.py
deleted file mode 100644
index bb6f558e5..000000000
--- a/src/calliope/backend/backends.py
+++ /dev/null
@@ -1,1367 +0,0 @@
-# Copyright (C) since 2013 Calliope contributors listed in AUTHORS.
-# Licensed under the Apache 2.0 License (see LICENSE file).
-
-from __future__ import annotations
-
-import logging
-import os
-import re
-import typing
-from abc import ABC, abstractmethod
-from contextlib import contextmanager, redirect_stderr, redirect_stdout
-from pathlib import Path
-from typing import (
-    Any,
-    Callable,
-    Generic,
-    Iterable,
-    Iterator,
-    Literal,
-    Optional,
-    TypeVar,
-    Union,
-)
-
-import numpy as np
-import pandas as pd
-import pyomo.environ as pe  # type: ignore
-import pyomo.kernel as pmo  # type: ignore
-import xarray as xr
-from pyomo.common.tempfiles import TempfileManager  # type: ignore
-from pyomo.opt import SolverFactory  # type: ignore
-
-from calliope.backend import parsing
-from calliope.core.util.logging import LogWriter
-from calliope.exceptions import BackendError, BackendWarning
-from calliope.exceptions import warn as model_warn
-
-T = TypeVar("T")
-_COMPONENTS_T = Literal[
-    "variables", "constraints", "objectives", "parameters", "global_expressions"
-]
-
-logger = logging.getLogger(__name__)
-
-
-class BackendModel(ABC, Generic[T]):
-    _VALID_COMPONENTS: tuple[_COMPONENTS_T, ...] = typing.get_args(_COMPONENTS_T)
-    _COMPONENT_ATTR_METADATA = ["description", "unit"]
-
-    def __init__(self, instance: T):
-        """Abstract base class for interfaces to solvers.
-
-        Args:
-            instance (T): Interface model instance.
-        """
-
-        self._instance = instance
-        self._dataset = xr.Dataset()
-        self.valid_math_element_names: set = set()
-
-    @abstractmethod
-    def add_parameter(
-        self,
-        parameter_name: str,
-        parameter_values: xr.DataArray,
-        default: Any = np.nan,
-        use_inf_as_na: bool = False,
-    ) -> None:
-        """
-        Add input parameter to backend model in-place.
-        If the backend interface allows for mutable parameter objects, they will be
-        generated, otherwise a copy of the model input dataset will be used.
-        In either case, NaN values are filled with the given parameter default value.
-
-        Args:
-            parameter_name (str): Name of parameter.
-            parameter_values (xr.DataArray): Array of parameter values.
-            default (Any, optional):
-                Default value to fill NaN entries in parameter values array.
-                Defaults to np.nan.
-            use_inf_as_na (bool, optional):
-                If True, will consider np.inf parameter value entries as np.nan and
-                consequently try to fill those entries with the parameter default value.
-                Defaults to False.
-        """
-
-    @abstractmethod
-    def add_constraint(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        constraint_dict: parsing.UnparsedConstraintDict,
-    ) -> None:
-        """
-        Add constraint equation to backend model in-place.
-        Resulting backend dataset entries will be constraint objects.
-
-        Args:
-            model_data (xr.Dataset):
-                Calliope model data with which to create an array mask - only those
-                dataset entries in the mask will be generated.
-            name (str):
-                Name of the constraint
-            constraint_dict (parsing.UnparsedConstraintDict):
-                Constraint configuration dictionary, ready to be parsed and then evaluated.
-        """
-
-    @abstractmethod
-    def add_global_expression(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        expression_dict: parsing.UnparsedExpressionDict,
-    ) -> None:
-        """
-        Add global expression (arithmetic combination of parameters and/or decision variables)
-        to backend model in-place.
-        Resulting backend dataset entries will be linear expression objects.
-
-        Args:
-            model_data (xr.Dataset):
-                Calliope model data with which to create an array mask - only those
-                dataset entries in the mask will be generated.
-            name (str):
-                Name of the global expression
-            expression_dict (parsing.UnparsedExpressionDict):
-                Global expression configuration dictionary, ready to be parsed and then evaluated.
-        """
-
-    @abstractmethod
-    def add_variable(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        variable_dict: parsing.UnparsedVariableDict,
-    ) -> None:
-        """
-        Add decision variable to backend model in-place.
-        Resulting backend dataset entries will be decision variable objects.
-
-        Args:
-            model_data (xr.Dataset):
-                Calliope model data with which to create an array mask - only those
-                dataset entries in the mask will be generated.
-            name (str):
-                Name of the variable.
-            variable_dict (parsing.UnparsedVariableDict):
-                Variable configuration dictionary, ready to be parsed and then evaluated.
-        """
-
-    @abstractmethod
-    def add_objective(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        objective_dict: parsing.UnparsedObjectiveDict,
-    ) -> None:
-        """
-        Add objective arithmetic to backend model in-place.
-        Resulting backend dataset entry will be a single, unindexed objective object.
-
-        Args:
-            model_data (xr.Dataset):
-                Calliope model data with which to create a constraint mask - only those
-                dataset entries in the mask will be generated.
-            name (str):
-                Name of the objective.
-            objective_dict (parsing.UnparsedObjectiveDict):
-                Objective configuration dictionary, ready to be parsed and then evaluated.
-        """
-
-    @abstractmethod
-    def get_parameter(
-        self, parameter_name: str, as_backend_objs: bool = True
-    ) -> Optional[xr.DataArray]:
-        """
-        Extract parameter from backend dataset.
-
-        Args:
-            parameter_name (str): Name of parameter.
-            TODO: hide this and create a method to edit parameter values
-                  (to handle interfaces with non-mutable params)
-            as_backend_objs (bool, optional):
-                If True, will keep the array entries as backend interface objects,
-                which can be updated to update the underlying model.
-                Otherwise, parameter values are given directly, with default values in place of NaNs.
-                Defaults to True.
-
-        Returns:
-            Optional[xr.DataArray]: If parameter is not in backend dataset, will return None.
-        """
-
-    @abstractmethod
-    def get_constraint(
-        self,
-        constraint_name: str,
-        as_backend_objs: bool = True,
-        eval_body: bool = False,
-    ) -> Optional[Union[xr.DataArray, xr.Dataset]]:
-        """
-        Get constraint data as either a table of details or as an array of backend interface objects.
-        Can be used to inspect and debug built constraints.
-
-        Args:
-            constraint_name (str): Name of constraint, as given in YAML constraint key.
-            TODO: hide this and create a method to edit constraints that handles differences in interface APIs.
-            as_backend_objs (bool, optional):
-                If True, will keep the array entries as backend interface objects,
-                which can be updated to change the underlying model.
-                Otherwise, constraint body, and lower and upper bounds are given in a table.
-                Defaults to True.
-            eval_body (bool, optional):
-                If True and as_backend_objs is False, will attempt to evaluate the constraint body.
-                If the model has been optimised, this attempt will produce a numeric value to see where the constraint sits between the lower or upper bound.
-                If the model has not yet been optimised, this attempt will fall back on the same as
-                if `eval_body` was set to False, i.e. a string representation of the linear expression in the constraint body.
-                Defaults to False.
-
-        Returns:
-            Optional[Union[xr.DataArray, xr.Dataset]]:
-                If constraint is not in backend dataset, will return None.
-                If as_backend_objs is True, will return an xr.DataArray.
-                Otherwise, a xr.Dataset will be given, indexed over the same dimensions as the xr.DataArray, with variables for the constraint body, and upper (`ub`) and lower (`lb`) bounds.
-        """
-
-    @abstractmethod
-    def get_variable(
-        self, variable_name: str, as_backend_objs: bool = True
-    ) -> Optional[xr.DataArray]:
-        """Extract decision variable array from backend dataset
-
-        Args:
-            variable_name (str): Name of variable.
-            TODO: hide this and create a method to edit variables that handles differences in interface APIs.
-            as_backend_objs (bool, optional):
-                If True, will keep the array entries as backend interface objects,
-                which can be updated to update the underlying model.
-                Otherwise, variable values are given directly.
-                If the model has not been successfully optimised, variable values will all be None.
-                Defaults to True.
-
-        Returns:
-            Optional[xr.DataArray]: If decision variable is not in backend dataset, will return None.
-        """
-
-    @abstractmethod
-    def get_global_expression(
-        self, expression_name: str, as_backend_objs: bool = True, eval_body: bool = True
-    ) -> Optional[xr.DataArray]:
-        """Extract global expression array from backend dataset
-
-        Args:
-            global_expression_name (str): Name of global expression
-            TODO: hide this and create a method to edit expressions that handles differences in interface APIs.
-            as_backend_objs (bool, optional):
-                If True, will keep the array entries as backend interface objects,
-                which can be updated to update the underlying model.
-                Otherwise, global expression values are given directly.
-                If the model has not been successfully optimised, expression values will all be provided as strings.
-                Defaults to True.
-            eval_body (bool, optional):
-                If True and as_backend_objs is False, will attempt to evaluate the expression.
-                If the model has been optimised, this attempt will produce a numeric value.
-                If the model has not yet been optimised, this attempt will fall back on the same as
-                if `eval_body` was set to False, i.e. a string representation of the linear expression.
-                Defaults to True.
-
-        Returns:
-            Optional[xr.DataArray]: If global_expression is not in backend dataset, will return None.
-        """
-
-    @abstractmethod
-    def solve(
-        self,
-        solver: str,
-        solver_io: Optional[str] = None,
-        solver_options: Optional[dict] = None,
-        save_logs: Optional[str] = None,
-        warmstart: bool = False,
-        **solve_kwargs,
-    ):
-        """
-        Optimise built model. If solution is optimal, interface objects
-        (decision variables, global expressions, constraints, objective) can be successfully
-        evaluated for their values at optimality.
-
-        Args:
-            solver (str): Name of solver to optimise with.
-            solver_io (Optional[str], optional):
-                If chosen solver has a python interface, set to "python" for potential
-                performance gains, otherwise should be left as None. Defaults to None.
-            solver_options (Optional[dict], optional):
-                Solver options/parameters to pass directly to solver.
-                See solver documentation for available parameters that can be influenced.
-                Defaults to None.
-            save_logs (Optional[str], optional):
-                If given, solver logs and built LP file will be saved to this filepath.
-                Defaults to None.
-            warmstart (bool, optional):
-                If True, and the chosen solver is capable of implementing it, an existing
-                optimal solution will be used to warmstart the next solve run.
-                Defaults to False.
-        """
-
-    @abstractmethod
-    def create_obj_list(self, key: str, component_type: _COMPONENTS_T) -> None:
-        """Attach an empty list object to the backend model object.
-        This may be a backend-specific subclass of a standard list object.
-
-        Args:
-            key (str): Name of object.
-            component_type (str): Object type.
-
-        Raises:
-            BackendError: Cannot overwrite object of same name and type.
-        """
-
-    @abstractmethod
-    def delete_obj_list(self, key: str, component_type: _COMPONENTS_T) -> None:
-        """Delete a list object from the backend model object.
-
-        Args:
-            key (str): Name of object.
-            component_type (str): Object type.
-        """
-
-    def load_results(self) -> xr.Dataset:
-        """
-        Evaluate backend decision variables, global expressions, and parameters (if not in inputs)
-        after a successful model run.
-
-        Returns:
-            xr.Dataset: Dataset of optimal solution results (all numeric data).
-        """
-
-        def _drop_attrs(da):
-            da.attrs = {
-                k: v for k, v in da.attrs.items() if k in self._COMPONENT_ATTR_METADATA
-            }
-            return da
-
-        all_variables = {
-            name_: _drop_attrs(self.get_variable(name_, as_backend_objs=False))
-            for name_, var in self.variables.items()
-            if var.notnull().any()
-        }
-        all_global_expressions = {
-            name_: _drop_attrs(
-                self.get_global_expression(name_, as_backend_objs=False, eval_body=True)
-            )
-            for name_, expr in self.global_expressions.items()
-            if expr.notnull().any()
-        }
-
-        results = xr.Dataset({**all_variables, **all_global_expressions})
-
-        return results
-
-    def add_all_parameters(self, model_data: xr.Dataset, run_config: dict) -> None:
-        """
-        Add all parameters to backend dataset in-place, including those in the run configuration.
-        If model data does not include a parameter, their default values will be added here
-        as unindexed backend dataset parameters.
-
-        TODO: Move the decision on which run config params to generate as backend params
-              earlier in the process.
-        Parameters in "objective_options" and the bigM parameter will be added from
-        run configuration.
-
-        Args:
-            model_data (xr.Dataset): Input model data.
-            defaults (dict): Parameter defaults.
-            run_config (UpdateObserverDict): Run configuration dictionary.
-        """
-
-        for param_name, param_data in model_data.filter_by_attrs(
-            is_result=0
-        ).data_vars.items():
-            default_val = model_data.attrs["defaults"].get(param_name, np.nan)
-            self.add_parameter(param_name, param_data, default_val)
-        for param_name, default_val in model_data.attrs["defaults"].items():
-            if param_name in self.parameters.keys():
-                continue
-            self.add_parameter(
-                param_name, xr.DataArray(default_val), use_inf_as_na=False
-            )
-
-        for option_name, option_val in run_config.get("objective_options", {}).items():
-            if option_name == "cost_class":
-                objective_cost_class = {
-                    k: v for k, v in option_val.items() if k in model_data.costs
-                }
-                self.add_parameter(
-                    "objective_cost_class",
-                    xr.DataArray.from_series(
-                        pd.Series(objective_cost_class).rename_axis(index="costs")
-                    ),
-                )
-            else:
-                self.add_parameter("objective_" + option_name, xr.DataArray(option_val))
-        self.add_parameter("bigM", xr.DataArray(run_config.get("bigM", 1e10)))
-
-    def apply_func(
-        self, func: Callable, *args, output_core_dims: tuple = ((),), **kwargs
-    ) -> xr.DataArray:
-        """
-        Apply a function to every element of an arbitrary number of xarray DataArrays.
-
-        Args:
-            func (Callable):
-                Un-vectorized function to call.
-                Number of accepted args should equal len(args).
-                Number of accepted kwargs should equal len(kwargs).
-            args (xr.DataArray):
-                xarray DataArrays which will be broadcast together and then iterated over
-                to apply the function.
-            output_core_dims (tuple):
-                Additional dimensions which are expected to be passed back from `xr.apply_ufunc` after applying `func`.
-                This is directly passed to `xr.apply_ufunc`; see their documentation for more details.
-                Defaults to ((), )
-            kwargs (dict[str, Any]):
-                Additional keyword arguments to pass to `func`.
-
-        Returns:
-            xr.DataArray: Array with func applied to all elements.
-        """
-        return xr.apply_ufunc(
-            func,
-            *args,
-            kwargs=kwargs,
-            vectorize=True,
-            keep_attrs=True,
-            output_dtypes=[np.dtype("O")],
-            output_core_dims=output_core_dims,
-        )
-
-    @abstractmethod
-    def verbose_strings(self) -> None:
-        """
-        Update optimisation model object string representations to include the index coordinates of the object.
-
-        E.g., `variables(carrier_prod)[0]` will become `variables(carrier_prod)[power, region1, ccgt, 2005-01-01 00:00]`
-
-        This takes approximately 10% of the peak memory required to initially build the optimisation problem, so should only be invoked if inspecting the model in detail (e.g., debugging)
-
-        Only string representations of model parameters and variables will be updated since global expressions automatically show the string representation of their contents.
-        """
-
-    @abstractmethod
-    def to_lp(self, path: Union[str, Path]) -> None:
-        """Write the optimisation problem to file in the linear programming LP format.
-        The LP file can be used for debugging and to submit to solvers directly.
-
-        Args:
-            path (Union[str, Path]): Path to which the LP file will be written.
-        """
-
-    def _raise_error_on_preexistence(self, key: str, obj_type: _COMPONENTS_T):
-        """
-        We do not allow any overlap of backend object names since they all have to
-        co-exist in the backend dataset.
-        I.e., users cannot overwrite any backend component with another
-        (of the same type or otherwise).
-
-        Args:
-            key (str): Backend object name
-            obj_type (Literal["variables", "constraints", "objectives", "parameters", "expressions"]): Object type.
-
-        Raises:
-            BackendError:
-                Raised if `key` already exists in the backend model
-                (either with the same or different type as `obj_type`).
-        """
-        if key in self._dataset.keys():
-            if key in getattr(self, obj_type):
-                raise BackendError(
-                    f"Trying to add already existing `{key}` to backend model {obj_type}."
-                )
-            else:
-                other_obj_type = self._dataset[key].attrs["obj_type"].removesuffix("s")
-                raise BackendError(
-                    f"Trying to add already existing *{other_obj_type}* `{key}` "
-                    f"as a backend model *{obj_type.removesuffix('s')}*."
-                )
-
-    @staticmethod
-    def _clean_arrays(*args) -> None:
-        """
-        Preemptively delete objects with large memory footprints that might otherwise
-        stick around longer than necessary.
-        """
-        del args
-
-    def _add_to_dataset(
-        self,
-        name: str,
-        da: xr.DataArray,
-        obj_type: _COMPONENTS_T,
-        unparsed_dict: Union[parsing.UNPARSED_DICTS, dict],
-        references: Optional[set] = None,
-    ):
-        """
-        Add array of backend objects to backend dataset in-place.
-
-        Args:
-            name (str): Name of entry in dataset.
-            da (xr.DataArray): Data to add.
-            obj_type (str): Type of backend objects in the array.
-            unparsed_dict (DT):
-                Dictionary describing the object being added, from which descriptor attributes will be extracted and added to the array attributes.
-            references (set):
-                All other backend objects which are references in this backend object's linear expression(s).
-                E.g. the constraint "carrier_prod / energy_eff <= energy_cap" references the variables ["carrier_prod", "energy_cap"]
-                and the parameter ["energy_eff"].
-                All referenced objects will have their "references" attribute updated with this object's name.
-                Defaults to None.
-        """
-
-        add_attrs = {
-            attr: unparsed_dict.get(attr)
-            for attr in self._COMPONENT_ATTR_METADATA
-            if attr in unparsed_dict.keys()
-        }
-
-        da.attrs.update(
-            {
-                "obj_type": obj_type,
-                "references": set(),
-                "coords_in_name": False,
-                **add_attrs,  # type: ignore
-            }
-        )
-        self._dataset[name] = da
-
-        if references is not None:
-            for reference in references:
-                self._dataset[reference].attrs["references"].add(name)
-
-    @property
-    def constraints(self):
-        "Slice of backend dataset to show only built constraints"
-        return self._dataset.filter_by_attrs(obj_type="constraints")
-
-    @property
-    def variables(self):
-        "Slice of backend dataset to show only built variables"
-        return self._dataset.filter_by_attrs(obj_type="variables")
-
-    @property
-    def parameters(self):
-        "Slice of backend dataset to show only built parameters"
-        return self._dataset.filter_by_attrs(obj_type="parameters")
-
-    @property
-    def global_expressions(self):
-        "Slice of backend dataset to show only built global expressions"
-        return self._dataset.filter_by_attrs(obj_type="global_expressions")
-
-    @property
-    def objectives(self):
-        "Slice of backend dataset to show only built objectives"
-        return self._dataset.filter_by_attrs(obj_type="objectives")
-
-
-class PyomoBackendModel(BackendModel):
-    _COMPONENT_TRANSLATOR = {
-        "parameter": "parameter",
-        "variable": "variable",
-        "global_expression": "expression",
-        "constraint": "constraint",
-        "objective": "objective",
-    }
-
-    def __init__(self, **kwargs):
-        BackendModel.__init__(self, instance=pmo.block())
-        self._instance.parameters = pmo.parameter_dict()
-        self._instance.variables = pmo.variable_dict()
-        self._instance.global_expressions = pmo.expression_dict()
-        self._instance.constraints = pmo.constraint_dict()
-        self._instance.objectives = pmo.objective_list()
-
-    def add_parameter(
-        self,
-        parameter_name: str,
-        parameter_values: xr.DataArray,
-        default: Any = np.nan,
-        use_inf_as_na: bool = False,
-    ) -> None:
-        self._raise_error_on_preexistence(parameter_name, "parameters")
-
-        self.create_obj_list(parameter_name, "parameters")
-
-        parameter_da = self.apply_func(
-            self._to_pyomo_param,
-            parameter_values,
-            name=parameter_name,
-            default=default,
-            use_inf_as_na=use_inf_as_na,
-        )
-        if parameter_da.isnull().all():
-            self.delete_obj_list(parameter_name, "parameters")
-            parameter_da = parameter_da.astype(float)
-
-        parameter_da.attrs["original_dtype"] = parameter_values.dtype
-        self._add_to_dataset(parameter_name, parameter_da, "parameters", {})
-        self.valid_math_element_names.add(parameter_name)
-
-    def add_constraint(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        constraint_dict: parsing.UnparsedConstraintDict,
-    ) -> None:
-        def _constraint_setter(
-            where: xr.DataArray, expr: tuple[xr.DataArray, str, xr.DataArray]
-        ) -> xr.DataArray:
-            lhs, op, rhs = expr
-            lhs = lhs.squeeze(drop=True)
-            rhs = rhs.squeeze(drop=True)
-
-            self._check_expr_where_consistency(lhs, where, f"(constraints, {name})")
-            self._check_expr_where_consistency(rhs, where, f"(constraints, {name})")
-
-            to_fill = self.apply_func(
-                self._to_pyomo_constraint,
-                where,
-                lhs,
-                rhs,
-                op=op,
-                name=name,
-            )
-            return to_fill
-
-        self._add_constraint_or_global_expression(
-            model_data,
-            name,
-            constraint_dict,
-            _constraint_setter,
-            "constraints",
-        )
-
-    def add_global_expression(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        expression_dict: parsing.UnparsedExpressionDict,
-    ) -> None:
-        def _expression_setter(where: xr.DataArray, expr: xr.DataArray) -> xr.DataArray:
-            expr = expr.squeeze(drop=True)
-
-            self._check_expr_where_consistency(expr, where, f"(expressions, {name})")
-
-            to_fill = self.apply_func(
-                self._to_pyomo_expression,
-                where,
-                expr,
-                name=name,
-            )
-            self._clean_arrays(expr)
-            return to_fill
-
-        self.valid_math_element_names.add(name)
-
-        self._add_constraint_or_global_expression(
-            model_data,
-            name,
-            expression_dict,
-            _expression_setter,
-            "global_expressions",
-        )
-
-    def add_variable(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        variable_dict: parsing.UnparsedVariableDict,
-    ) -> None:
-        self.valid_math_element_names.add(name)
-
-        parsed_variable = parsing.ParsedBackendComponent(
-            "variables", name, variable_dict
-        )
-
-        where = parsed_variable.generate_top_level_where_array(model_data)
-        if not where.any():
-            return None
-
-        self._raise_error_on_preexistence(name, "variables")
-        self.create_obj_list(name, "variables")
-
-        domain = parsed_variable._unparsed.get("domain", "real")
-        domain_type = getattr(pmo, f"{domain.title()}Set")
-
-        lb, ub = self._get_capacity_bounds(variable_dict["bounds"], name=name)
-        variable_da = self.apply_func(
-            self._to_pyomo_variable,
-            where,
-            ub,
-            lb,
-            name=name,
-            domain_type=domain_type,
-        )
-
-        self._add_to_dataset(name, variable_da, "variables", variable_dict)
-
-    def add_objective(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        objective_dict: parsing.UnparsedObjectiveDict,
-    ) -> None:
-        self._raise_error_on_preexistence(name, "objectives")
-        sense_dict = {"minimize": 1, "minimise": 1, "maximize": -1, "maximise": -1}
-        parsed_objective = parsing.ParsedBackendComponent(
-            "objectives", name, objective_dict
-        )
-        equations = parsed_objective.parse_equations(self.valid_math_element_names)
-
-        n_valid_exprs = 0
-        for equation in equations:
-            where = equation.evaluate_where(model_data)
-            if where.any():
-                expr = equation.evaluate_expression(
-                    model_data, self, where=where
-                ).item()
-                n_valid_exprs += 1
-
-        if n_valid_exprs == 0:
-            return None
-
-        if n_valid_exprs > 1:
-            raise BackendError(
-                f"More than one {name} objective is valid for this "
-                "optimisation problem; only one is allowed."
-            )
-
-        objective = pmo.objective(expr, sense=sense_dict[objective_dict["sense"]])
-
-        if name == model_data.run_config["objective"]:
-            objective.activate()
-        else:
-            objective.deactivate()
-
-        self._instance.objectives.append(objective)
-
-        self._add_to_dataset(
-            name, xr.DataArray(objective), "objectives", objective_dict
-        )
-
-    def get_parameter(
-        self, parameter_name: str, as_backend_objs: bool = True
-    ) -> Optional[xr.DataArray]:
-        parameter = self.parameters.get(parameter_name, None)
-        if as_backend_objs or not isinstance(parameter, xr.DataArray):
-            return parameter
-
-        param_as_vals = self.apply_func(self._from_pyomo_param, parameter)
-        if parameter.original_dtype.kind == "M":  # i.e., np.datetime64
-            return xr.apply_ufunc(pd.to_datetime, param_as_vals)
-        else:
-            return param_as_vals.astype(parameter.original_dtype)
-
-    def get_constraint(
-        self,
-        constraint_name: str,
-        as_backend_objs: bool = True,
-        eval_body: bool = False,
-    ) -> Optional[Union[xr.DataArray, xr.Dataset]]:
-        constraint = self.constraints.get(constraint_name, None)
-        if isinstance(constraint, xr.DataArray) and not as_backend_objs:
-            constraint_attrs = self.apply_func(
-                self._from_pyomo_constraint,
-                constraint,
-                eval_body=eval_body,
-                output_core_dims=(["attributes"],),
-            )
-            constraint_attrs.coords["attributes"] = ["lb", "body", "ub"]
-            constraint = constraint_attrs.to_dataset("attributes")
-        return constraint
-
-    def get_variable(
-        self, variable_name: str, as_backend_objs: bool = True
-    ) -> Optional[xr.DataArray]:
-        variable = self.variables.get(variable_name, None)
-        if isinstance(variable, xr.DataArray) and not as_backend_objs:
-            return self.apply_func(self._from_pyomo_param, variable)
-        else:
-            return variable
-
-    def get_global_expression(
-        self,
-        global_expression_name: str,
-        as_backend_objs: bool = True,
-        eval_body: bool = False,
-    ) -> Optional[xr.DataArray]:
-        global_expression = self.global_expressions.get(global_expression_name, None)
-        if isinstance(global_expression, xr.DataArray) and not as_backend_objs:
-            return self.apply_func(
-                self._from_pyomo_expr, global_expression, eval_body=eval_body
-            )
-        else:
-            return global_expression
-
-    def solve(
-        self,
-        solver: str,
-        solver_io: Optional[str] = None,
-        solver_options: Optional[dict] = None,
-        save_logs: Optional[str] = None,
-        warmstart: bool = False,
-        **solve_kwargs,
-    ):
-        opt = SolverFactory(solver, solver_io=solver_io)
-
-        if solver_options:
-            for k, v in solver_options.items():
-                opt.options[k] = v
-
-        if save_logs is not None:
-            solve_kwargs.update({"symbolic_solver_labels": True, "keepfiles": True})
-            os.makedirs(save_logs, exist_ok=True)
-            TempfileManager.tempdir = save_logs  # Sets log output dir
-        if warmstart and solver in ["glpk", "cbc"]:
-            model_warn(
-                "The chosen solver, {}, does not support warmstart, which may "
-                "impact performance.".format(solver)
-            )
-            warmstart = False
-
-        with redirect_stdout(LogWriter(logger, "debug", strip=True)):  # type: ignore
-            with redirect_stderr(LogWriter(logger, "error", strip=True)):  # type: ignore
-                # Ignore most of gurobipy's logging, as it's output is
-                # already captured through STDOUT
-                logging.getLogger("gurobipy").setLevel(logging.ERROR)
-                results = opt.solve(self._instance, tee=True, **solve_kwargs)
-
-        termination = results.solver[0].termination_condition
-
-        if termination == pe.TerminationCondition.optimal:
-            self._instance.load_solution(results.solution[0])
-
-        else:
-            logger.critical("Problem status:")
-            for line in str(results.problem[0]).split("\n"):
-                logger.critical(line)
-            logger.critical("Solver status:")
-            for line in str(results.solver[0]).split("\n"):
-                logger.critical(line)
-
-            model_warn("Model solution was non-optimal.", _class=BackendWarning)
-
-        return str(termination)
-
-    def verbose_strings(self) -> None:
-        def __renamer(val, *idx):
-            if pd.notnull(val):
-                val.calliope_coords = idx
-
-        with self._datetime_as_string(self._dataset):
-            for component_type in ["parameters", "variables"]:
-                for da in self._dataset.filter_by_attrs(
-                    coords_in_name=False, **{"obj_type": component_type}
-                ).values():
-                    self.apply_func(__renamer, da, *[da.coords[i] for i in da.dims])
-                    da.attrs["coords_in_name"] = True
-
-    def to_lp(self, path: Union[str, Path]) -> None:
-        self._instance.write(str(path), format="lp", symbolic_solver_labels=True)
-
-    def create_obj_list(self, key: str, component_type: _COMPONENTS_T) -> None:
-        """Attach an empty pyomo kernel list object to the pyomo model object.
-
-        Args:
-            key (str): Name of object
-            component_type (str): Object type
-
-        Raises:
-            BackendError: Cannot overwrite object of same name and type.
-        """
-        component_dict = getattr(self._instance, component_type)
-        if key in component_dict:
-            raise BackendError(
-                f"Trying to add already existing `{key}` to backend model {component_type}."
-            )
-        else:
-            singular_component = component_type.removesuffix("s")
-            component_dict[key] = getattr(
-                pmo, f"{self._COMPONENT_TRANSLATOR[singular_component]}_list"
-            )()
-
-    def delete_obj_list(self, key: str, component_type: _COMPONENTS_T) -> None:
-        """Delete a list object from the backend model object.
-
-        Args:
-            key (str): Name of object
-            component_type (str): Object type
-        """
-        component_dict = getattr(self._instance, component_type)
-        if key not in component_dict:
-            return None
-        else:
-            del component_dict[key]
-
-    def _add_constraint_or_global_expression(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        component_dict: Union[
-            parsing.UnparsedConstraintDict, parsing.UnparsedExpressionDict
-        ],
-        component_setter: Callable,
-        component_type: Literal["constraints", "global_expressions"],
-    ) -> None:
-        """Generalised function to add a constraint or global expression array to the model.
-
-        Args:
-            model_data (xr.Dataset): Calliope model input data
-            name: Name of the constraint or global expression
-            component_dict (Union[parsing.UnparsedConstraintDict, parsing.UnparsedExpressionDict]):
-                Unparsed YAML dictionary configuration.
-            component_setter (Callable):
-                Function to combine evaluated xarray DataArrays into
-                constraint/global expression objects.
-                Will receive outputs of `evaluate_where` and `evaluate_expression` as inputs.
-            component_type (Literal[constraints, global_expressions])
-            parser (Callable): Parsing rule to use for the component (differs between constraints and global expressions)
-
-
-        Raises:
-            BackendError:
-                The sub-equations of the parsed component cannot generate component
-                objects on duplicate index entries.
-        """
-        references: set[str] = set()
-        parsed_component = parsing.ParsedBackendComponent(
-            component_type, name, component_dict
-        )
-
-        top_level_where = parsed_component.generate_top_level_where_array(
-            model_data, align_to_foreach_sets=False
-        )
-        if not top_level_where.any():
-            return None
-
-        self._raise_error_on_preexistence(name, component_type)
-        component_da = (
-            xr.DataArray()
-            .where(parsed_component.drop_dims_not_in_foreach(top_level_where))
-            .astype(np.dtype("O"))
-        )
-        self.create_obj_list(name, component_type)
-
-        equations = parsed_component.parse_equations(self.valid_math_element_names)
-        for element in equations:
-            where = element.evaluate_where(model_data, initial_where=top_level_where)
-            if not where.any():
-                continue
-
-            where = parsed_component.drop_dims_not_in_foreach(where)
-
-            if component_da.where(where).notnull().any():
-                subset_overlap = component_da.where(where).to_series().dropna().index
-
-                raise BackendError(
-                    "Trying to set two equations for the same index of "
-                    f"{component_type.removesuffix('s')} `{name}`:\n{subset_overlap}"
-                )
-
-            expr = element.evaluate_expression(
-                model_data, self, where=where, references=references
-            )
-            to_fill = component_setter(where, expr)
-            component_da = component_da.fillna(to_fill)
-
-        if component_da.isnull().all():
-            self.delete_obj_list(name, component_type)
-            return None
-
-        self._add_to_dataset(
-            name, component_da, component_type, component_dict, references
-        )
-
-    @staticmethod
-    def _check_expr_where_consistency(
-        expression: xr.DataArray, where: xr.DataArray, description: str
-    ) -> None:
-        """
-        Checks if a given constraint or global expression is consistent with the binary where array.
-
-        Parameters:
-            expression (xr.DataArray): array of linear expressions from a global expression or one side of a constraint equation.
-            where (xr.DataArray): where array.
-            description (str): Description to prefix the error message.
-
-        Raises:
-            BackendError:
-                Raised if there is a dimension in the expression that is not in the where array.
-            BackendError:
-                Raised if the expression has any NaN where the where array applies.
-        """
-        # Check whether expression has a dim that does not exist in where.
-        broadcast_dims_where = set(expression.dims).difference(set(where.dims))
-        if broadcast_dims_where:
-            raise BackendError(
-                f"{description}: The linear expression array is indexed over dimensions not present in `foreach`: {broadcast_dims_where}"
-            )
-
-        incomplete_constraints = expression.isnull() & where
-        if incomplete_constraints.any():
-            raise BackendError(
-                f"{description}: Missing a linear expression for some coordinates selected by 'where'. Adapting 'where' might help."
-            )
-
-    def _get_capacity_bounds(
-        self, bounds: parsing.UnparsedVariableBoundDict, name: str
-    ) -> tuple[xr.DataArray, xr.DataArray]:
-        """
-        Generate arrays corresponding to upper and lower bounds of a decision variable.
-        If `equals` is given, then everywhere it is not None/np.nan it will be applied
-        as the simultaneous upper and lower bound. Everywhere it is None/np.nan, it will
-        be filled by `min` (for lower bound) and `max` (for upper bound).
-        Upper and lower bounds will be scaled by `scale`, if `scale` is not None/np.nan.
-
-        Args:
-            bounds (dict): Dictionary of optional keys `min`, `max`, `equals`, and `scale`.
-            name (str): Name of decision variable.
-
-        Returns:
-            tuple[xr.DataArray, xr.DataArray]:
-                (upper bounds, lower bounds). Where unbounded, the array entry will be None.
-        """
-
-        def __get_bound(bound):
-            this_bound = bounds.get(bound, None)
-            if isinstance(this_bound, str):
-                return self.get_parameter(this_bound)
-            else:
-                # TODO: decide if this parameter should be added to the backend dataset too
-                name_ = f"TEMP_{name}_{bound}"
-                self.create_obj_list(name_, "parameters")
-                return xr.DataArray(self._to_pyomo_param(this_bound, name=name_))
-
-        scale = __get_bound("scale")
-        equals_ = __get_bound("equals")
-        min_ = __get_bound("min")
-        max_ = __get_bound("max")
-
-        lb = equals_.fillna(min_)
-        ub = equals_.fillna(max_)
-        if scale.notnull().any():
-            lb = lb * scale
-            ub = ub * scale
-
-        return lb.fillna(None), ub.fillna(None)
-
-    def _to_pyomo_param(
-        self, val: Any, *, name: str, default: Any = np.nan, use_inf_as_na: bool = True
-    ) -> Union[type[ObjParameter], float]:
-        """
-        Utility function to generate a pyomo parameter for every element of an
-        xarray DataArray.
-        Output objects are of the type ObjParameter(pmo.parameter) since they need a
-        "dtype" property to be handled by xarray.
-
-        If not np.nan/None, output objects are also added to the backend model object in-place.
-
-
-        Args:
-            val (Any): Value to turn into a mutable pyomo parameter
-            name (str): Name of parameter
-            default (Any, optional): Default value if `val` is None/np.nan. Defaults to np.nan.
-            use_inf_as_na (bool, optional): If True, see np.inf as np.nan. Defaults to True.
-
-        Returns:
-            Union[type[ObjParameter], float]:
-                If both `val` and `default` are np.nan/None, return np.nan.
-                Otherwise return ObjParameter(val/default).
-        """
-        if use_inf_as_na:
-            val = np.nan if val in [np.inf, -np.inf] else val
-            default = np.nan if default in [np.inf, -np.inf] else default
-
-        if pd.isnull(val):
-            if pd.isnull(default):
-                param = np.nan
-            else:
-                param = ObjParameter(default)
-                self._instance.parameters[name].append(param)
-        else:
-            param = ObjParameter(val)
-            self._instance.parameters[name].append(param)
-        return param
-
-    def _to_pyomo_constraint(
-        self,
-        mask: Union[bool, np.bool_],
-        lhs: Any,
-        rhs: Any,
-        *,
-        op: Literal["==", ">=", "<="],
-        name: str,
-    ) -> Union[type[pmo.constraint], float]:
-        """
-        Utility function to generate a pyomo constraint for every element of an
-        xarray DataArray.
-
-        If not np.nan/None, output objects are also added to the backend model object in-place.
-
-        Args:
-            mask (Union[bool, np.bool_]): If True, add constraint, otherwise return np.nan
-            lhs (Any): Equation left-hand-side linear expression
-            rhs (Any): Equation right-hand-side linear expression
-
-        Kwargs:
-            op (Literal[, optional): Operator to compare `lhs` and `rhs`. Defaults to =", ">=", "<="].
-            name (str): Name of constraint
-
-        Returns:
-            Union[type[pmo.constraint], float]:
-                If mask is True, return np.nan.
-                Otherwise return pmo_constraint(expr=lhs op rhs).
-        """
-
-        if not mask:
-            return np.nan
-        elif op == "==":
-            constraint = pmo.constraint(expr=lhs == rhs)
-        elif op == "<=":
-            constraint = pmo.constraint(expr=lhs <= rhs)
-        elif op == ">=":
-            constraint = pmo.constraint(expr=lhs >= rhs)
-        self._instance.constraints[name].append(constraint)
-        return constraint
-
-    def _to_pyomo_expression(
-        self, mask: Union[bool, np.bool_], expr: Any, *, name: str
-    ) -> Union[type[pmo.expression], float]:
-        """
-        Utility function to generate a pyomo expression for every element of an
-        xarray DataArray.
-
-        If not np.nan/None, output objects are also added to the backend model object in-place.
-
-
-        Args:
-            mask (Union[bool, np.bool_]): If True, add expression, otherwise return np.nan.
-            expr (Any): Linear expression to add.
-        Kwargs:
-            name (str): Expression name.
-
-        Returns:
-            Union[type[pmo.expression], float]:
-                If mask is True, return np.nan.
-                Otherwise return pmo_expression(expr).
-        """
-        if mask:
-            expr_obj = pmo.expression(expr)
-            self._instance.global_expressions[name].append(expr_obj)
-            return expr_obj
-        else:
-            return np.nan
-
-    def _to_pyomo_variable(
-        self,
-        mask: Union[bool, np.bool_],
-        ub: Any,
-        lb: Any,
-        *,
-        name: str,
-        domain_type: Literal["RealSet", "IntegerSet"],
-    ) -> Union[type[ObjVariable], float]:
-        """
-        Utility function to generate a pyomo decision variable for every element of an
-        xarray DataArray.
-
-        If not np.nan/None, output objects are also added to the backend model object in-place.
-
-        Args:
-            mask (Union[bool, np.bool_]): If True, add variable, otherwise return np.nan.
-            ub (Any): Upper bound to apply to the variable.
-            lb (Any): Lower bound to apply to the variable.
-
-        Kwargs:
-            domain_type (Literal["RealSet", "IntegerSet"]):
-                Domain over which variables are valid (real = continuous, integer = integer/binary)
-            name (str): Name of variable.
-
-        Returns:
-            Union[type[ObjVariable], float]:
-                If mask is True, return np.nan.
-                Otherwise return pmo_variable(ub=ub, lb=lb, domain_type=domain_type).
-        """
-        if mask:
-            var = ObjVariable(ub=ub, lb=lb, domain_type=domain_type)
-            self._instance.variables[name].append(var)
-            return var
-        else:
-            return np.nan
-
-    @staticmethod
-    def _from_pyomo_param(val: Union[ObjParameter, ObjVariable, float]) -> Any:
-        """
-        Evaluate value of Pyomo object.
-        If the input object is a parameter, a numeric/string value will be given.
-        If the input object is a global expression or variable, a numeric value will be given
-        only if the backend model has been successfully optimised, otherwise evaluation will return None.
-
-        Args:
-            val (Union[ObjParameter, pmo.expression, ObjVariable, np.nan]):
-                Item to be evaluated.
-
-        Returns:
-            Any: If the input is nullable, return np.nan, otherwise evaluate the "value" method of the object.
-        """
-        if pd.isnull(val):
-            return np.nan
-        else:
-            return val.value  # type: ignore
-
-    @staticmethod
-    def _from_pyomo_constraint(
-        val: pmo.constraint, *, eval_body: bool = False
-    ) -> pd.Series:
-        """Evaluate Pyomo constraint object.
-
-        Args:
-            val (pmo.constraint): constraint object to be evaluated
-        Kwargs:
-            eval_body (bool, optional):
-                If True, attempt to evaluate the constraint object `body`, which will evaluate the
-                linear expression contained in the constraint body and produce a numeric value.
-                This will only succeed if the backend model has been successfully optimised,
-                otherwise a string representation of the linear expression will be returned
-                (same as eval_body=False). Defaults to False.
-
-        Returns:
-            pd.Series: Array of upper bound (ub), body, and lower bound (lb).
-        """
-        if pd.isnull(val):
-            vals = [np.nan, np.nan, np.nan]
-        else:
-            if eval_body:
-                try:
-                    body = val.body()
-                except ValueError:
-                    body = val.body.to_string()
-            else:
-                body = val.body.to_string()
-            vals = [val.lb, body, val.ub]
-        return pd.Series(data=vals, index=["lb", "body", "ub"])
-
-    @staticmethod
-    def _from_pyomo_expr(val: pmo.expression, *, eval_body: bool = False) -> Any:
-        """Evaluate Pyomo expression object.
-
-        Args:
-            val (pmo.expression): expression object to be evaluated
-            eval_body (bool, optional):
-                If True, attempt to evaluate the expression object, which will produce a numeric value.
-                This will only succeed if the backend model has been successfully optimised,
-                otherwise a string representation of the linear expression will be returned
-                (same as eval_body=False). Defaults to False.
-
-        Returns:
-            Any: If the input is nullable, return np.nan, otherwise a numeric value
-            (eval_body=True and problem is optimised) or a string.
-        """
-        if pd.isnull(val):
-            return np.nan
-        else:
-            if eval_body:
-                expr = val()
-                if expr is None:
-                    return val.to_string()
-                else:
-                    return expr
-            else:
-                return val.to_string()
-
-    @contextmanager
-    def _datetime_as_string(self, data: Union[xr.DataArray, xr.Dataset]) -> Iterator:
-        """Context manager to temporarily convert np.dtype("datetime64[ns]") coordinates (e.g. timesteps) to strings with a resolution of minutes.
-
-        Args:
-            data (Union[xr.DataArray, xr.Dataset]): xarray object on whose coordinates the conversion will take place.
-        """
-        datetime_coords = set()
-        for name_, vals_ in data.coords.items():
-            if vals_.dtype.kind == "M":
-                data.coords[name_] = data.coords[name_].dt.strftime("%Y-%m-%d %H:%M")
-                datetime_coords.add(name_)
-        try:
-            yield
-        finally:
-            for name_ in datetime_coords:
-                data.coords[name_] = xr.apply_ufunc(
-                    pd.to_datetime, data.coords[name_], keep_attrs=True
-                )
-
-
-class CoordObj(ABC):
-    """Class with methods to update the `name` property of inheriting classes"""
-
-    def __init__(self) -> None:
-        self._calliope_coords: Optional[Iterable] = None
-
-    def _update_name(self, old_name: str) -> str:
-        """
-        Update string of a list containing a single number with a string of a list containing any arbitrary number of elements
-
-        Args:
-            old_name (str): String representation of a list containing a single number
-
-        Returns:
-            str:
-                If `self.calliope_coords` is None, returns `old_name`.
-                Otherwise returns string representation of a list containing the contents of `self.calliope_coords`
-        """
-        if self._calliope_coords is None:
-            return old_name
-
-        if not self._calliope_coords:  # empty list = dimensionless component
-            coord_list = ""
-        else:
-            coord_list = f"[{', '.join(str(i) for i in self._calliope_coords)}]"
-        return re.sub(r"\[\d+\]", coord_list, old_name)
-
-    @property
-    def calliope_coords(self):
-        return self._calliope_coords
-
-    @calliope_coords.setter
-    def calliope_coords(self, val):
-        self._calliope_coords = val
-
-
-class ObjParameter(pmo.parameter, CoordObj):
-    """
-    A pyomo parameter (`a object for storing a mutable, numeric value that can be used to build a symbolic expression`)
-    with added `dtype` property and a `name` property setter (via the `pmo.parameter.getname` method) which replaces a list position as a name with a list of strings.
-    """
-
-    def __init__(self, value, **kwds):
-        assert not pd.isnull(value)
-        pmo.parameter.__init__(self, value, **kwds)
-        CoordObj.__init__(self)
-
-    @property
-    def dtype(self):
-        return "O"
-
-    def getname(self, *args, **kwargs):
-        return self._update_name(pmo.parameter.getname(self, *args, **kwargs))
-
-
-class ObjVariable(pmo.variable, CoordObj):
-    """
-    A pyomo variable with a `name` property setter (via the `pmo.variable.getname` method) which replaces a list position as a name with a list of strings.
-
-    """
-
-    def __init__(self, **kwds):
-        pmo.variable.__init__(self, **kwds)
-        CoordObj.__init__(self)
-
-    def getname(self, *args, **kwargs):
-        return self._update_name(pmo.variable.getname(self, *args, **kwargs))
diff --git a/src/calliope/backend/expression_parser.py b/src/calliope/backend/expression_parser.py
index 4417f2bcd..cf0bee6a7 100644
--- a/src/calliope/backend/expression_parser.py
+++ b/src/calliope/backend/expression_parser.py
@@ -51,7 +51,7 @@
 from calliope.exceptions import BackendError
 
 if TYPE_CHECKING:
-    from calliope.backend.backends import BackendModel
+    from calliope.backend.backend_model import BackendModel
 
 pp.ParserElement.enablePackrat()
 
diff --git a/src/calliope/backend/latex_backend.py b/src/calliope/backend/latex_backend_model.py
similarity index 61%
rename from src/calliope/backend/latex_backend.py
rename to src/calliope/backend/latex_backend_model.py
index 89fbc2331..d86a40617 100644
--- a/src/calliope/backend/latex_backend.py
+++ b/src/calliope/backend/latex_backend_model.py
@@ -9,8 +9,8 @@
 import numpy as np
 import xarray as xr
 
-from calliope.backend import backends, parsing
-from calliope.exceptions import BackendError, ModelError
+from calliope.backend import backend_model, parsing
+from calliope.exceptions import ModelError
 
 _ALLOWED_MATH_FILE_FORMATS = Literal["tex", "rst", "md"]
 
@@ -21,9 +21,10 @@ def __init__(self, backend_builder: Callable) -> None:
 
         Args:
             backend_builder (Callable):
-                Method to generate all optimisation problem components on a calliope.backends.BackendModel object.
+                Method to generate all optimisation problem components on a calliope.backend_model.BackendModel object.
         """
         self._builder = backend_builder
+        self._inputs: xr.Dataset
 
     def build(
         self,
@@ -36,10 +37,18 @@ def build(
                 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".
         """
 
-        backend = LatexBackendModel(include=include)
+        backend = LatexBackendModel(self._inputs, include=include)
 
         self._instance = self._builder(backend)
 
+    @property
+    def inputs(self):
+        return self._inputs
+
+    @inputs.setter
+    def inputs(self, val: xr.Dataset):
+        self._inputs = val
+
     @overload  # noqa: F811
     def write(  # noqa: F811
         self,
@@ -101,30 +110,30 @@ def write(  # noqa: F811
             return None
 
 
-class LatexBackendModel(backends.BackendModel):
+class LatexBackendModel(backend_model.BackendModelGenerator):
     # \negthickspace used to counter the introduction of spaces to separate the curly braces
     # in \text. Curly braces need separating otherwise jinja2 gets confused.
     LATEX_EQUATION_ELEMENT = textwrap.dedent(
         r"""
         \begin{array}{r}
-        {% if sets is defined and sets%}
+        {% if sets is defined and sets %}
             \forall{}
         {% for set in sets %}
             \text{ {{set|removesuffix("s")}} }\negthickspace \in \negthickspace\text{ {{set + "," if not loop.last else set }} }
         {% endfor %}
             \\
         {% endif %}
-        {% if sense is defined %}
+        {% if sense is defined and sense %}
             {{sense}}
         {% endif %}
-        {% if where is defined and where != "" %}
+        {% if where is defined and where and where != "" %}
             \text{if } {{where}}
         {% endif %}
         \end{array}
         \begin{cases}
         {% for equation in equations %}
             {{equation["expression"]}}&\quad
-        {% if "where" in equation and equation["where"] != "" %}
+        {% if "where" in equation and equation.where != "" %}
             \text{if } {{equation["where"]}}
         {% endif %}
             \\
@@ -246,16 +255,20 @@ class LatexBackendModel(backends.BackendModel):
     )
     FORMAT_STRINGS = {"rst": RST_DOC, "tex": TEX_DOC, "md": MD_DOC}
 
-    def __init__(self, include: Literal["all", "valid"] = "all") -> None:
+    def __init__(
+        self, inputs: xr.Dataset, include: Literal["all", "valid"] = "all"
+    ) -> None:
         """Interface to build a string representation of the mathematical formulation using LaTeX math notation.
 
         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".
         """
-        backends.BackendModel.__init__(self, instance=dict())
+        super().__init__(inputs)
         self.include = include
 
+        self._add_all_inputs_as_parameters()
+
     def add_parameter(
         self,
         parameter_name: str,
@@ -268,148 +281,129 @@ def add_parameter(
 
     def add_constraint(
         self,
-        model_data: xr.Dataset,
         name: str,
-        constraint_dict: parsing.UnparsedConstraintDict,
+        constraint_dict: Optional[parsing.UnparsedConstraintDict] = None,
     ) -> None:
-        self._add_constraint_or_expression(
-            model_data,
+        equation_strings: list = []
+
+        def _constraint_setter(
+            element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set
+        ) -> xr.DataArray:
+            self._add_latex_strings(where, element, equation_strings)
+            return where.where(where)
+
+        parsed_component = self._add_component(
             name,
             constraint_dict,
-            lambda x: None,
+            _constraint_setter,
             "constraints",
+            break_early=False,
+        )
+
+        self._generate_math_string(
+            parsed_component,
+            self.constraints[name],
+            equations=equation_strings,
         )
 
     def add_global_expression(
         self,
-        model_data: xr.Dataset,
         name: str,
-        expression_dict: parsing.UnparsedConstraintDict,
+        expression_dict: Optional[parsing.UnparsedExpressionDict] = None,
     ) -> None:
         self.valid_math_element_names.add(name)
 
-        self._add_constraint_or_expression(
-            model_data,
+        equation_strings: list = []
+
+        def _expression_setter(
+            element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set
+        ) -> xr.DataArray:
+            self._add_latex_strings(where, element, equation_strings)
+            return where.where(where)
+
+        parsed_component = self._add_component(
             name,
             expression_dict,
-            lambda x: None,
+            _expression_setter,
             "global_expressions",
+            break_early=False,
+        )
+
+        self._generate_math_string(
+            parsed_component,
+            self.global_expressions[name],
+            equations=equation_strings,
         )
 
     def add_variable(
         self,
-        model_data: xr.Dataset,
         name: str,
-        variable_dict: parsing.UnparsedVariableDict,
+        variable_dict: Optional[parsing.UnparsedVariableDict] = None,
     ) -> None:
+        domain_dict = {"real": r"\R", "integer": r"\Z"}
+
         self.valid_math_element_names.add(name)
-        self._raise_error_on_preexistence(name, "variables")
 
-        parsed_variable = parsing.ParsedBackendComponent(
-            "variables", name, variable_dict
-        )
-        where_array = parsed_variable.generate_top_level_where_array(
-            model_data, break_early=False
-        )
+        if variable_dict is None:
+            variable_dict = self.inputs.math["variables"][name]
 
-        # add early to be accessed when creating bound strings.
-        self._add_to_dataset(name, where_array, "variables", variable_dict)
+        parsed_component = self._add_component(
+            name, variable_dict, lambda where: where, "variables", break_early=False
+        )
+        where_array = self.variables[name]
 
-        where_latex = parsed_variable.evaluate_where(model_data, as_latex=True)
-        lb, ub = self._get_capacity_bounds(variable_dict["bounds"], name, model_data)
+        domain = domain_dict[variable_dict.get("domain", "real")]
+        lb, ub = self._get_capacity_bounds(name, variable_dict["bounds"])
 
-        if self.include == "all" or (self.include == "valid" and where_array.any()):
-            self._generate_math_string(
-                where_array,
-                sets=parsed_variable.sets,
-                where=where_latex if where_latex != "" else None,
-                equations=[lb, ub],
-            )
-        # add again to ensure "math_string" attribute is there.
-        self._add_to_dataset(name, where_array, "variables", variable_dict)
+        self._generate_math_string(
+            parsed_component,
+            where_array,
+            equations=[lb, ub],
+            sense=r"\forall" + domain,
+        )
 
     def add_objective(
         self,
-        model_data: xr.Dataset,
         name: str,
-        objective_dict: parsing.UnparsedObjectiveDict,
+        objective_dict: Optional[parsing.UnparsedObjectiveDict] = None,
     ) -> None:
-        self._raise_error_on_preexistence(name, "objectives")
         sense_dict = {
             "minimize": r"\min{}",
             "maximize": r"\max{}",
             "minimise": r"\min{}",
             "maximise": r"\max{}",
         }
-        parsed_objective = parsing.ParsedBackendComponent(
-            "objectives", name, objective_dict
-        )
-        equations = parsed_objective.parse_equations(self.valid_math_element_names)
-        equation_strings = []
-        for element in equations:
-            if self.include == "valid":
-                where_array = element.evaluate_where(model_data)
-            where_latex = element.evaluate_where(model_data, as_latex=True)
-            expr = element.evaluate_expression(model_data, self, as_latex=True)
-
-            if self.include == "all" or (self.include == "valid" and where_array.any()):
-                equation_strings.append({"expression": expr, "where": where_latex})
-        objective_da = xr.DataArray()
-        if equation_strings:
-            self._generate_math_string(
-                objective_da,
-                sense=sense_dict[objective_dict["sense"]],
-                equations=equation_strings,
-            )
-        self._add_to_dataset(name, objective_da, "objectives", objective_dict)
-
-    def get_parameter(
-        self, parameter_name: str, as_backend_objs: bool = True
-    ) -> Optional[xr.DataArray]:
-        return self.parameters.get(parameter_name, None)
-
-    def create_obj_list(self, key: str, component_type: backends._COMPONENTS_T) -> None:
-        return None
+        if objective_dict is None:
+            objective_dict = self.inputs.math["objectives"][name]
+        equation_strings: list = []
+
+        def _objective_setter(
+            element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set
+        ) -> None:
+            self._add_latex_strings(where, element, equation_strings)
+            return None
 
-    def delete_obj_list(self, key: str, component_type: backends._COMPONENTS_T) -> None:
-        return None
+        parsed_component = self._add_component(
+            name, objective_dict, _objective_setter, "objectives", break_early=False
+        )
 
-    def get_constraint(
-        self,
-        constraint_name: str,
-        as_backend_objs: bool = True,
-        eval_body: bool = False,
-    ) -> Optional[Union[xr.DataArray, xr.Dataset]]:
-        return self.constraints.get(constraint_name, None)
-
-    def get_variable(
-        self, variable_name: str, as_backend_objs: bool = True
-    ) -> Optional[xr.DataArray]:
-        return self.variables.get(variable_name, None)
-
-    def get_global_expression(
-        self, expression_name: str, as_backend_objs: bool = True, eval_body: bool = True
-    ) -> Optional[xr.DataArray]:
-        return self.global_expressions.get(expression_name, None)
-
-    def solve(
-        self,
-        solver: str,
-        solver_io: Optional[str] = None,
-        solver_options: Optional[dict] = None,
-        save_logs: Optional[str] = None,
-        warmstart: bool = False,
-        **solve_kwargs,
-    ):
-        raise BackendError(
-            "Cannot solve a LaTex backend model - this only exists to produce a string representation of the model math"
+        self._generate_math_string(
+            parsed_component,
+            self.objectives[name],
+            equations=equation_strings,
+            sense=sense_dict[objective_dict["sense"]],
         )
 
-    def verbose_strings(self):
+    def _create_obj_list(
+        self, key: str, component_type: backend_model._COMPONENTS_T
+    ) -> None:
         return None
 
-    def to_lp(self, path: Union[str, Path]):
-        return None
+    def delete_component(
+        self, key: str, component_type: backend_model._COMPONENTS_T
+    ) -> None:
+        if key in self._dataset and self._dataset[key].obj_type == component_type:
+            del self._dataset[key]
 
     def generate_math_doc(self, format: _ALLOWED_MATH_FILE_FORMATS = "tex") -> str:
         """Generate the math documentation by embedding LaTeX math in a template.
@@ -442,64 +436,37 @@ def generate_math_doc(self, format: _ALLOWED_MATH_FILE_FORMATS = "tex") -> str:
         }
         return self._render(doc_template, components=components)
 
-    def _add_constraint_or_expression(
-        self,
-        model_data: xr.Dataset,
-        name: str,
-        component_dict: Union[
-            parsing.UnparsedConstraintDict, parsing.UnparsedExpressionDict
-        ],
-        component_setter: Callable,
-        component_type: Literal["constraints", "global_expressions"],
-    ) -> None:
-        parsed_component = parsing.ParsedBackendComponent(
-            component_type, name, component_dict
-        )
+    def _add_latex_strings(self, where, element, equation_strings):
+        expr = element.evaluate_expression(self.inputs, self, as_latex=True)
+        where_latex = element.evaluate_where(self.inputs, as_latex=True)
 
-        top_level_where = parsed_component.generate_top_level_where_array(
-            model_data, break_early=False
-        )
-        component_da = xr.DataArray().where(
-            parsed_component.drop_dims_not_in_foreach(top_level_where)
-        )
-        top_level_where_latex = parsed_component.evaluate_where(
-            model_data, as_latex=True
-        )
-
-        self._raise_error_on_preexistence(name, component_type)
-
-        equations = parsed_component.parse_equations(self.valid_math_element_names)
-        equation_strings = []
-        for element in equations:
-            where_array = element.evaluate_where(
-                model_data, initial_where=top_level_where
-            )
-            where_array = parsed_component.drop_dims_not_in_foreach(where_array)
-
-            expr = element.evaluate_expression(model_data, self, as_latex=True)
-            where_latex = element.evaluate_where(model_data, as_latex=True)
-
-            if self.include == "all" or (self.include == "valid" and where_array.any()):
-                equation_strings.append({"expression": expr, "where": where_latex})
-            component_da = component_da.fillna(where_array.where(where_array))
-        if equation_strings:
-            self._generate_math_string(
-                component_da,
-                sets=parsed_component.sets,
-                where=top_level_where_latex,
-                equations=equation_strings,
-            )
-        self._add_to_dataset(
-            name, component_da.fillna(0), component_type, component_dict
-        )
+        if self.include == "all" or (self.include == "valid" and where.any()):
+            equation_strings.append({"expression": expr, "where": where_latex})
 
     def _generate_math_string(
         self,
-        da: xr.DataArray,
-        **kwargs,
+        parsed_component: Optional[parsing.ParsedBackendComponent],
+        where_array: xr.DataArray,
+        equations: Optional[list[dict[str, str]]] = None,
+        sense: Optional[str] = None,
+        where: Optional[str] = None,
+        sets: Optional[list[str]] = None,
     ) -> None:
-        equation_element_string = self._render(self.LATEX_EQUATION_ELEMENT, **kwargs)
-        da.attrs.update({"math_string": equation_element_string})
+        if parsed_component is not None:
+            where = parsed_component.evaluate_where(self.inputs, as_latex=True)
+            sets = parsed_component.sets
+
+        if self.include == "all" or (
+            self.include == "valid" and where_array.fillna(0).any()
+        ):
+            equation_element_string = self._render(
+                self.LATEX_EQUATION_ELEMENT,
+                equations=equations if equations is not None else [],
+                sense=sense,
+                where=where,
+                sets=sets,
+            )
+            where_array.attrs.update({"math_string": equation_element_string})
         return None
 
     @staticmethod
@@ -510,13 +477,10 @@ def _render(template: str, **kwargs) -> str:
 
     def _get_capacity_bounds(
         self,
-        bounds: parsing.UnparsedVariableBoundDict,
         name: str,
-        model_data: xr.Dataset,
+        bounds: parsing.UnparsedVariableBoundDict,
     ) -> tuple[dict[str, str], ...]:
         bound_dict: parsing.UnparsedConstraintDict = {
-            "foreach": [],
-            "where": "True",
             "equations": [
                 {"expression": f"{bounds['min']} <= {name}"},
                 {"expression": f"{name} <= {bounds['max']}"},
@@ -527,6 +491,6 @@ def _get_capacity_bounds(
             self.valid_math_element_names,
         )
         return tuple(
-            {"expression": eq.evaluate_expression(model_data, self, as_latex=True)}
+            {"expression": eq.evaluate_expression(self.inputs, self, as_latex=True)}
             for eq in equations
         )
diff --git a/src/calliope/backend/parsing.py b/src/calliope/backend/parsing.py
index c49f27d3e..6362ecc64 100644
--- a/src/calliope/backend/parsing.py
+++ b/src/calliope/backend/parsing.py
@@ -6,14 +6,28 @@
 import functools
 import itertools
 import operator
-from typing import Any, Callable, Iterable, Literal, Optional, TypeVar, Union, overload
+from typing import (
+    TYPE_CHECKING,
+    Any,
+    Callable,
+    Iterable,
+    Literal,
+    Optional,
+    TypeVar,
+    Union,
+    overload,
+)
 
 import pyparsing as pp
 import xarray as xr
 from typing_extensions import NotRequired, Required, TypedDict
 
 from calliope import exceptions
-from calliope.backend import backends, expression_parser, helper_functions, where_parser
+from calliope.backend import expression_parser, helper_functions, where_parser
+
+if TYPE_CHECKING:
+    from calliope.backend import backend_model
+
 
 TRUE_ARRAY = xr.DataArray(True)
 
@@ -25,8 +39,8 @@ class UnparsedEquationDict(TypedDict):
 
 class UnparsedConstraintDict(TypedDict):
     description: NotRequired[str]
-    foreach: Required[list]
-    where: str
+    foreach: NotRequired[list]
+    where: NotRequired[str]
     equations: Required[list[UnparsedEquationDict]]
     sub_expressions: NotRequired[dict[str, list[UnparsedEquationDict]]]
     slices: NotRequired[dict[str, list[UnparsedEquationDict]]]
@@ -296,7 +310,7 @@ def drop_dims_not_in_foreach(self, where: xr.DataArray) -> xr.DataArray:
     def evaluate_expression(  # noqa: F811
         self,
         model_data: xr.Dataset,
-        backend_interface: backends.BackendModel,
+        backend_interface: backend_model.BackendModel,
         as_latex: Literal[False] = False,
         references: Optional[set] = None,
         where: Optional[xr.DataArray] = None,
@@ -307,7 +321,7 @@ def evaluate_expression(  # noqa: F811
     def evaluate_expression(  # noqa: F811
         self,
         model_data: xr.Dataset,
-        backend_interface: backends.BackendModel,
+        backend_interface: backend_model.BackendModel,
         as_latex: Literal[True],
         references: Optional[set] = None,
     ) -> str:
@@ -316,7 +330,7 @@ def evaluate_expression(  # noqa: F811
     def evaluate_expression(  # noqa: F811
         self,
         model_data: xr.Dataset,
-        backend_interface: backends.BackendModel,
+        backend_interface: backend_model.BackendModel,
         as_latex: bool = False,
         references: Optional[set] = None,
         where: Optional[xr.DataArray] = None,
diff --git a/src/calliope/backend/pyomo_backend_model.py b/src/calliope/backend/pyomo_backend_model.py
new file mode 100644
index 000000000..7be47211d
--- /dev/null
+++ b/src/calliope/backend/pyomo_backend_model.py
@@ -0,0 +1,914 @@
+# Copyright (C) since 2013 Calliope contributors listed in AUTHORS.
+# Licensed under the Apache 2.0 License (see LICENSE file).
+
+from __future__ import annotations
+
+import logging
+import os
+import re
+from abc import ABC
+from contextlib import contextmanager, redirect_stderr, redirect_stdout
+from pathlib import Path
+from typing import (
+    Any,
+    Iterable,
+    Iterator,
+    Literal,
+    Optional,
+    SupportsFloat,
+    TypeVar,
+    Union,
+)
+
+import numpy as np
+import pandas as pd
+import pyomo.environ as pe  # type: ignore
+import pyomo.kernel as pmo  # type: ignore
+import xarray as xr
+from pyomo.common.tempfiles import TempfileManager  # type: ignore
+from pyomo.opt import SolverFactory  # type: ignore
+
+from calliope.backend import backend_model, parsing
+from calliope.core.util.logging import LogWriter
+from calliope.exceptions import BackendError, BackendWarning
+from calliope.exceptions import warn as model_warn
+
+T = TypeVar("T")
+_COMPONENTS_T = Literal[
+    "variables", "constraints", "objectives", "parameters", "global_expressions"
+]
+
+logger = logging.getLogger("calliope.backend.backend_model")
+
+COMPONENT_TRANSLATOR = {
+    "parameter": "parameter",
+    "variable": "variable",
+    "global_expression": "expression",
+    "constraint": "constraint",
+    "objective": "objective",
+}
+
+
+class PyomoBackendModel(backend_model.BackendModel):
+    def __init__(self, inputs: xr.Dataset):
+        super().__init__(inputs, instance=pmo.block())
+
+        self._instance.parameters = pmo.parameter_dict()
+        self._instance.variables = pmo.variable_dict()
+        self._instance.global_expressions = pmo.expression_dict()
+        self._instance.constraints = pmo.constraint_dict()
+        self._instance.objectives = pmo.objective_dict()
+
+        self._add_all_inputs_as_parameters()
+
+    def add_parameter(
+        self,
+        parameter_name: str,
+        parameter_values: xr.DataArray,
+        default: Any = np.nan,
+        use_inf_as_na: bool = False,
+    ) -> None:
+        self._raise_error_on_preexistence(parameter_name, "parameters")
+
+        self._create_obj_list(parameter_name, "parameters")
+
+        parameter_da = self._apply_func(
+            self._to_pyomo_param,
+            parameter_values,
+            name=parameter_name,
+            default=default,
+            use_inf_as_na=use_inf_as_na,
+        )
+        if parameter_da.isnull().all():
+            self.delete_component(parameter_name, "parameters")
+            parameter_da = parameter_da.astype(float)
+
+        parameter_da.attrs["original_dtype"] = parameter_values.dtype
+        self._add_to_dataset(parameter_name, parameter_da, "parameters", {})
+        self.valid_math_element_names.add(parameter_name)
+
+    def add_constraint(
+        self,
+        name: str,
+        constraint_dict: Optional[parsing.UnparsedConstraintDict] = None,
+    ) -> None:
+        def _constraint_setter(
+            element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set
+        ) -> xr.DataArray:
+            lhs, op, rhs = element.evaluate_expression(
+                self.inputs, self, where=where, references=references
+            )
+            lhs = lhs.squeeze(drop=True)
+            rhs = rhs.squeeze(drop=True)
+
+            self._check_expr_where_consistency(lhs, where, f"(constraints, {name})")
+            self._check_expr_where_consistency(rhs, where, f"(constraints, {name})")
+
+            to_fill = self._apply_func(
+                self._to_pyomo_constraint,
+                where,
+                lhs,
+                rhs,
+                op=op,
+                name=name,
+            )
+            return to_fill
+
+        self._add_component(name, constraint_dict, _constraint_setter, "constraints")
+
+    def add_global_expression(
+        self,
+        name: str,
+        expression_dict: Optional[parsing.UnparsedExpressionDict] = None,
+    ) -> None:
+        def _expression_setter(
+            element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set
+        ) -> xr.DataArray:
+            expr = element.evaluate_expression(
+                self.inputs, self, where=where, references=references
+            )
+            expr = expr.squeeze(drop=True)
+
+            self._check_expr_where_consistency(expr, where, f"(expressions, {name})")
+
+            to_fill = self._apply_func(
+                self._to_pyomo_expression,
+                where,
+                expr,
+                name=name,
+            )
+            self._clean_arrays(expr)
+            return to_fill
+
+        self.valid_math_element_names.add(name)
+
+        self._add_component(
+            name, expression_dict, _expression_setter, "global_expressions"
+        )
+
+    def add_variable(
+        self,
+        name: str,
+        variable_dict: Optional[parsing.UnparsedVariableDict] = None,
+    ) -> None:
+        domain_dict = {"real": pmo.RealSet, "integer": pmo.IntegerSet}
+
+        if variable_dict is None:
+            variable_dict = self.inputs.math["variables"][name]
+        self.valid_math_element_names.add(name)
+
+        def _variable_setter(where):
+            domain_type = domain_dict[variable_dict.get("domain", "real")]
+
+            lb, ub = self._get_capacity_bounds(name, variable_dict["bounds"])
+            return self._apply_func(
+                self._to_pyomo_variable,
+                where,
+                ub,
+                lb,
+                name=name,
+                domain_type=domain_type,
+            )
+
+        self.valid_math_element_names.add(name)
+
+        self._add_component(name, variable_dict, _variable_setter, "variables")
+
+    def add_objective(
+        self,
+        name: str,
+        objective_dict: Optional[parsing.UnparsedObjectiveDict] = None,
+    ) -> None:
+        sense_dict = {"minimize": 1, "minimise": 1, "maximize": -1, "maximise": -1}
+
+        if objective_dict is None:
+            objective_dict = self.inputs.math["objectives"][name]
+        sense = sense_dict[objective_dict["sense"]]
+
+        def _objective_setter(
+            element: parsing.ParsedBackendEquation, where: xr.DataArray, references: set
+        ) -> xr.DataArray:
+            expr = element.evaluate_expression(self.inputs, self, references=references)
+            objective = pmo.objective(expr.item(), sense=sense)
+            if name == self.inputs.run_config["objective"]:
+                objective.activate()
+            else:
+                objective.deactivate()
+
+            self._instance.objectives[name].append(objective)
+            return xr.DataArray(objective)
+
+        self._add_component(name, objective_dict, _objective_setter, "objectives")
+
+    def get_parameter(
+        self,
+        name: str,
+        as_backend_objs: bool = True,
+    ) -> xr.DataArray:
+        parameter = self.parameters.get(name, None)
+        if parameter is None:
+            raise KeyError(f"Unknown parameter: {name}")
+
+        if as_backend_objs or not isinstance(parameter, xr.DataArray):
+            return parameter
+
+        param_as_vals = self._apply_func(self._from_pyomo_param, parameter)
+        if parameter.original_dtype.kind == "M":  # i.e., np.datetime64
+            return xr.apply_ufunc(pd.to_datetime, param_as_vals)
+        else:
+            return param_as_vals.astype(parameter.original_dtype)
+
+    def get_constraint(
+        self,
+        name: str,
+        as_backend_objs: bool = True,
+        eval_body: bool = False,
+    ) -> Union[xr.DataArray, xr.Dataset]:
+        constraint = self.constraints.get(name, None)
+        if constraint is None:
+            raise KeyError(f"Unknown constraint: {name}")
+        if isinstance(constraint, xr.DataArray) and not as_backend_objs:
+            constraint_attrs = self._apply_func(
+                self._from_pyomo_constraint,
+                constraint,
+                eval_body=eval_body,
+                output_core_dims=(["attributes"],),
+            )
+            constraint_attrs.coords["attributes"] = ["lb", "body", "ub"]
+            constraint = constraint_attrs.to_dataset("attributes")
+        return constraint
+
+    def get_variable(
+        self,
+        name: str,
+        as_backend_objs: bool = True,
+    ) -> xr.DataArray:
+        variable = self.variables.get(name, None)
+        if variable is None:
+            raise KeyError(f"Unknown variable: {name}")
+        if as_backend_objs:
+            return variable
+        else:
+            return self._apply_func(self._from_pyomo_param, variable)
+
+    def get_variable_bounds(self, name: str) -> xr.Dataset:
+        variable = self.get_variable(name, as_backend_objs=True)
+        variable_attrs = self._apply_func(
+            self._from_pyomo_variable_bounds,
+            variable,
+            output_core_dims=(["attributes"],),
+        )
+        variable_attrs.coords["attributes"] = ["lb", "ub"]
+        return variable_attrs.to_dataset("attributes")
+
+    def get_global_expression(
+        self,
+        name: str,
+        as_backend_objs: bool = True,
+        eval_body: bool = False,
+    ) -> xr.DataArray:
+        global_expression = self.global_expressions.get(name, None)
+        if global_expression is None:
+            raise KeyError(f"Unknown global_expression: {name}")
+        if isinstance(global_expression, xr.DataArray) and not as_backend_objs:
+            return self._apply_func(
+                self._from_pyomo_expr, global_expression, eval_body=eval_body
+            )
+        else:
+            return global_expression
+
+    def solve(
+        self,
+        solver: str,
+        solver_io: Optional[str] = None,
+        solver_options: Optional[dict] = None,
+        save_logs: Optional[str] = None,
+        warmstart: bool = False,
+        **solve_kwargs,
+    ):
+        opt = SolverFactory(solver, solver_io=solver_io)
+
+        if solver_options:
+            for k, v in solver_options.items():
+                opt.options[k] = v
+
+        if save_logs is not None:
+            solve_kwargs.update({"symbolic_solver_labels": True, "keepfiles": True})
+            os.makedirs(save_logs, exist_ok=True)
+            TempfileManager.tempdir = save_logs  # Sets log output dir
+        if warmstart and solver in ["glpk", "cbc"]:
+            model_warn(
+                "The chosen solver, {}, does not support warmstart, which may "
+                "impact performance.".format(solver)
+            )
+            warmstart = False
+
+        with redirect_stdout(LogWriter(logger, "debug", strip=True)):  # type: ignore
+            with redirect_stderr(LogWriter(logger, "error", strip=True)):  # type: ignore
+                # Ignore most of gurobipy's logging, as it's output is
+                # already captured through STDOUT
+                logging.getLogger("gurobipy").setLevel(logging.ERROR)
+                results = opt.solve(self._instance, tee=True, **solve_kwargs)
+
+        termination = results.solver[0].termination_condition
+
+        if termination == pe.TerminationCondition.optimal:
+            self._instance.load_solution(results.solution[0])
+
+        else:
+            logger.critical("Problem status:")
+            for line in str(results.problem[0]).split("\n"):
+                logger.critical(line)
+            logger.critical("Solver status:")
+            for line in str(results.solver[0]).split("\n"):
+                logger.critical(line)
+
+            model_warn("Model solution was non-optimal.", _class=BackendWarning)
+
+        return str(termination)
+
+    def verbose_strings(self) -> None:
+        def __renamer(val, *idx):
+            if pd.notnull(val):
+                val.calliope_coords = idx
+
+        with self._datetime_as_string(self._dataset):
+            for component_type in ["parameters", "variables"]:
+                for da in self._dataset.filter_by_attrs(
+                    coords_in_name=False, **{"obj_type": component_type}
+                ).values():
+                    self._apply_func(__renamer, da, *[da.coords[i] for i in da.dims])
+                    da.attrs["coords_in_name"] = True
+
+    def to_lp(self, path: Union[str, Path]) -> None:
+        self._instance.write(str(path), format="lp", symbolic_solver_labels=True)
+
+    def _create_obj_list(self, key: str, component_type: _COMPONENTS_T) -> None:
+        """Attach an empty pyomo kernel list object to the pyomo model object.
+
+        Args:
+            key (str): Name of object
+            component_type (str): Object type
+
+        Raises:
+            BackendError: Cannot overwrite object of same name and type.
+        """
+        component_dict = getattr(self._instance, component_type)
+        if key in component_dict:
+            raise BackendError(
+                f"Trying to add already existing `{key}` to backend model {component_type}."
+            )
+        else:
+            singular_component = component_type.removesuffix("s")
+            component_dict[key] = getattr(
+                pmo, f"{COMPONENT_TRANSLATOR[singular_component]}_list"
+            )()
+
+    def delete_component(self, key: str, component_type: _COMPONENTS_T) -> None:
+        """Delete a list object from the backend model object.
+
+        Args:
+            key (str): Name of object
+            component_type (str): Object type
+        """
+        component_dict = getattr(self._instance, component_type)
+        if key in component_dict:
+            del component_dict[key]
+
+        if key in self._dataset and self._dataset[key].obj_type == component_type:
+            del self._dataset[key]
+
+    def update_parameter(
+        self, name: str, new_values: Union[xr.DataArray, SupportsFloat]
+    ) -> None:
+        new_values = xr.DataArray(new_values)
+        parameter_da = self.get_parameter(name)
+        missing_dims_in_new_vals = set(parameter_da.dims).difference(new_values.dims)
+        missing_dims_in_orig_vals = set(new_values.dims).difference(parameter_da.dims)
+        refs_to_update: set = set()
+
+        if (
+            (not parameter_da.shape and new_values.shape)
+            or missing_dims_in_orig_vals
+            or (parameter_da.isnull() & new_values.notnull()).any()
+        ):
+            refs_to_update = self._find_all_references(parameter_da.attrs["references"])
+            if refs_to_update:
+                logger.warning(
+                    f"Defining values for a previously fully/partially undefined parameter. The optimisation problem components {refs_to_update} will be re-built."
+                )
+            self.delete_component(name, "parameters")
+            self.add_parameter(
+                name, new_values, default=self.inputs.defaults.get(name, np.nan)
+            )
+            for ref in refs_to_update:
+                self._rebuild_reference(ref)
+            return None
+
+        if missing_dims_in_new_vals:
+            logger.warning(
+                f"New values will be broadcast along the {missing_dims_in_new_vals} dimension(s)."
+            )
+
+        self._apply_func(self._update_pyomo_param, parameter_da, new_values)
+
+    def update_variable_bounds(
+        self,
+        name: str,
+        *,
+        min: Optional[Union[xr.DataArray, SupportsFloat]] = None,
+        max: Optional[Union[xr.DataArray, SupportsFloat]] = None,
+    ) -> None:
+        translator = {"min": "lb", "max": "ub"}
+        variable_da = self.get_variable(name)
+
+        for bound_name, new_bounds in {"min": min, "max": max}.items():
+            if new_bounds is None:
+                continue
+
+            existing_bound_param = self.inputs.math.get_key(
+                f"variables.{name}.bounds.{bound_name}", None
+            )
+            if existing_bound_param in self.parameters:
+                raise BackendError(
+                    f"Cannot update variable bounds that have been set by parameters. use `update_parameter('{existing_bound_param}')` to update the {bound_name} bound of {name}."
+                )
+
+            bound_da = xr.DataArray(new_bounds)
+            missing_dims_in_new_vals = set(variable_da.dims).difference(bound_da.dims)
+            if missing_dims_in_new_vals:
+                logger.warning(
+                    f"New `{bound_name}` bounds for variable `{name}` will be broadcast along the {missing_dims_in_new_vals} dimension(s)."
+                )
+            self._apply_func(
+                self._update_pyomo_variable,
+                variable_da,
+                xr.DataArray(new_bounds),
+                bound=translator[bound_name],
+            )
+
+    def fix_variable(self, name: str, where: Optional[xr.DataArray] = None) -> None:
+        variable_da = self.get_variable(name)
+        if where is not None:
+            variable_da = variable_da.where(where.fillna(0))
+        self._apply_func(self._fix_pyomo_variable, variable_da)
+
+    def unfix_variable(self, name: str, where: Optional[xr.DataArray] = None) -> None:
+        variable_da = self.get_variable(name)
+        if where is not None:
+            variable_da = variable_da.where(where.fillna(0))
+        self._apply_func(self._unfix_pyomo_variable, variable_da)
+
+    @staticmethod
+    def _check_expr_where_consistency(
+        expression: xr.DataArray, where: xr.DataArray, description: str
+    ) -> None:
+        """
+        Checks if a given constraint or global expression is consistent with the where.
+
+        Parameters:
+            expression (xr.DataArray): array of linear expressions from a global expression or one side of a constraint equation.
+            where (xr.DataArray): where array.
+            description (str): Description to prefix the error message.
+
+        Raises:
+            BackendError:
+                Raised if there is a dimension in the expression that is not in the where.
+            BackendError:
+                Raised if the expression has any NaN where the where applies.
+        """
+        # Check whether expression has a dim that does not exist in where.
+        broadcast_dims_where = set(expression.dims).difference(set(where.dims))
+        if broadcast_dims_where:
+            raise BackendError(
+                f"{description}: The linear expression array is indexed over dimensions not present in `foreach`: {broadcast_dims_where}"
+            )
+
+        incomplete_constraints = expression.isnull() & where
+        if incomplete_constraints.any():
+            raise BackendError(
+                f"{description}: Missing a linear expression for some coordinates selected by 'where'. Adapting 'where' might help."
+            )
+
+    def _get_capacity_bounds(
+        self, name: str, bounds: parsing.UnparsedVariableBoundDict
+    ) -> tuple[xr.DataArray, xr.DataArray]:
+        """
+        Generate arrays corresponding to upper and lower bounds of a decision variable.
+        If `equals` is given, then everywhere it is not None/np.nan it will be applied
+        as the simultaneous upper and lower bound. Everywhere it is None/np.nan, it will
+        be filled by `min` (for lower bound) and `max` (for upper bound).
+        Upper and lower bounds will be scaled by `scale`, if `scale` is not None/np.nan.
+
+        Args:
+            bounds (dict): Dictionary of optional keys `min`, `max`, `equals`, and `scale`.
+            name (str): Name of decision variable.
+
+        Returns:
+            tuple[xr.DataArray, xr.DataArray]:
+                (upper bounds, lower bounds). Where unbounded, the array entry will be None.
+        """
+
+        def __get_bound(bound):
+            this_bound = bounds.get(bound, None)
+            if isinstance(this_bound, str):
+                return self.get_parameter(this_bound)
+            else:
+                # TODO: decide if this parameter should be added to the backend dataset too
+                name_ = f"__{name}_{bound}"
+                self._create_obj_list(name_, "parameters")
+                return xr.DataArray(self._to_pyomo_param(this_bound, name=name_))
+
+        scale = __get_bound("scale")
+        equals_ = __get_bound("equals")
+        min_ = __get_bound("min")
+        max_ = __get_bound("max")
+
+        lb = equals_.fillna(min_)
+        ub = equals_.fillna(max_)
+        if scale.notnull().any():
+            lb = lb * scale
+            ub = ub * scale
+
+        return lb.fillna(None), ub.fillna(None)
+
+    def _to_pyomo_param(
+        self, val: Any, *, name: str, default: Any = np.nan, use_inf_as_na: bool = True
+    ) -> Union[type[ObjParameter], float]:
+        """
+        Utility function to generate a pyomo parameter for every element of an
+        xarray DataArray.
+        Output objects are of the type ObjParameter(pmo.parameter) since they need a
+        "dtype" property to be handled by xarray.
+
+        If not np.nan/None, output objects are also added to the backend model object in-place.
+
+
+        Args:
+            val (Any): Value to turn into a mutable pyomo parameter
+            name (str): Name of parameter
+            default (Any, optional): Default value if `val` is None/np.nan. Defaults to np.nan.
+            use_inf_as_na (bool, optional): If True, see np.inf as np.nan. Defaults to True.
+
+        Returns:
+            Union[type[ObjParameter], float]:
+                If both `val` and `default` are np.nan/None, return np.nan.
+                Otherwise return ObjParameter(val/default).
+        """
+        if use_inf_as_na:
+            val = np.nan if val in [np.inf, -np.inf] else val
+            default = np.nan if default in [np.inf, -np.inf] else default
+        if pd.isnull(val):
+            if pd.isnull(default):
+                param = np.nan
+            else:
+                param = ObjParameter(default)
+                self._instance.parameters[name].append(param)
+        else:
+            param = ObjParameter(val)
+            self._instance.parameters[name].append(param)
+        return param
+
+    def _update_pyomo_param(self, orig: ObjParameter, new: Any) -> None:
+        """Utility function to update pyomo parameter values in-place.
+
+        Args:
+            orig (ObjParameter): Pyomo parameter to update.
+            new (Any): Value with which to update the parameter.
+        """
+        if pd.notnull(new):
+            orig.value = new
+
+    def _update_pyomo_variable(
+        self, orig: ObjVariable, new: Any, *, bound: Literal["lb", "ub"]
+    ) -> None:
+        """Utility function to update pyomo variable bounds in-place.
+
+        Args:
+            orig (ObjVariable): Pyomo variable to update.
+        Keyword Args:
+            lb (Any): Value with which to update the lower bound of the variable.
+            ub (Any): Value with which to update the upper bound of the variable.
+        """
+        if pd.notnull(orig) and pd.notnull(new):
+            setattr(orig, bound, new)
+
+    def _fix_pyomo_variable(self, orig: ObjVariable) -> None:
+        """Utility function to fix a pyomo variable to its value in the optimisation model solution.
+        Fixed variables will be considered as parameters in the subsequent solve.
+
+        Args:
+            orig (ObjVariable): Pyomo variable to fix.
+
+        Raises:
+            BackendError: Can only fix variables if they have values assigned to them from an optimal solution.
+        """
+        if pd.isnull(orig):
+            return None
+        elif orig.value is None:
+            raise BackendError(
+                "Cannot fix variable values without already having solved the model successfully."
+            )
+        else:
+            orig.fix()
+
+    def _unfix_pyomo_variable(self, orig: ObjVariable) -> None:
+        """Utility function to unfix a pyomo variable so it can be considered a decision variable on the next solve.
+
+        Args:
+            orig (ObjVariable): Pyomo variable to unfix.
+
+        """
+        if pd.notnull(orig):
+            orig.unfix()
+
+    def _to_pyomo_constraint(
+        self,
+        mask: Union[bool, np.bool_],
+        lhs: Any,
+        rhs: Any,
+        *,
+        op: Literal["==", ">=", "<="],
+        name: str,
+    ) -> Union[type[pmo.constraint], float]:
+        """
+        Utility function to generate a pyomo constraint for every element of an
+        xarray DataArray.
+
+        If not np.nan/None, output objects are also added to the backend model object in-place.
+
+        Args:
+            mask (Union[bool, np.bool_]): If True, add constraint, otherwise return np.nan
+            lhs (Any): Equation left-hand-side linear expression
+            rhs (Any): Equation right-hand-side linear expression
+
+        Kwargs:
+            op (Literal[, optional): Operator to compare `lhs` and `rhs`. Defaults to =", ">=", "<="].
+            name (str): Name of constraint
+
+        Returns:
+            Union[type[pmo.constraint], float]:
+                If mask is True, return np.nan.
+                Otherwise return pmo_constraint(expr=lhs op rhs).
+        """
+
+        if not mask:
+            return np.nan
+        elif op == "==":
+            constraint = pmo.constraint(expr=lhs == rhs)
+        elif op == "<=":
+            constraint = pmo.constraint(expr=lhs <= rhs)
+        elif op == ">=":
+            constraint = pmo.constraint(expr=lhs >= rhs)
+        self._instance.constraints[name].append(constraint)
+        return constraint
+
+    def _to_pyomo_expression(
+        self, mask: Union[bool, np.bool_], expr: Any, *, name: str
+    ) -> Union[type[pmo.expression], float]:
+        """
+        Utility function to generate a pyomo expression for every element of an
+        xarray DataArray.
+
+        If not np.nan/None, output objects are also added to the backend model object in-place.
+
+
+        Args:
+            mask (Union[bool, np.bool_]): If True, add expression, otherwise return np.nan.
+            expr (Any): Linear expression to add.
+        Kwargs:
+            name (str): Expression name.
+
+        Returns:
+            Union[type[pmo.expression], float]:
+                If mask is True, return np.nan.
+                Otherwise return pmo_expression(expr).
+        """
+        if mask:
+            expr_obj = pmo.expression(expr)
+            self._instance.global_expressions[name].append(expr_obj)
+            return expr_obj
+        else:
+            return np.nan
+
+    def _to_pyomo_variable(
+        self,
+        mask: Union[bool, np.bool_],
+        ub: Any,
+        lb: Any,
+        *,
+        name: str,
+        domain_type: Literal["RealSet", "IntegerSet"],
+    ) -> Union[type[ObjVariable], float]:
+        """
+        Utility function to generate a pyomo decision variable for every element of an
+        xarray DataArray.
+
+        If not np.nan/None, output objects are also added to the backend model object in-place.
+
+        Args:
+            mask (Union[bool, np.bool_]): If True, add variable, otherwise return np.nan.
+            ub (Any): Upper bound to apply to the variable.
+            lb (Any): Lower bound to apply to the variable.
+
+        Kwargs:
+            domain_type (Literal["RealSet", "IntegerSet"]):
+                Domain over which variables are valid (real = continuous, integer = integer/binary)
+            name (str): Name of variable.
+
+        Returns:
+            Union[type[ObjVariable], float]:
+                If mask is True, return np.nan.
+                Otherwise return pmo_variable(ub=ub, lb=lb, domain_type=domain_type).
+        """
+        if mask:
+            var = ObjVariable(ub=ub, lb=lb, domain_type=domain_type)
+            self._instance.variables[name].append(var)
+            return var
+        else:
+            return np.nan
+
+    @staticmethod
+    def _from_pyomo_param(val: Union[ObjParameter, ObjVariable, float]) -> Any:
+        """
+        Evaluate value of Pyomo object.
+        If the input object is a parameter, a numeric/string value will be given.
+        If the input object is a global expression or variable, a numeric value will be given
+        only if the backend model has been successfully optimised, otherwise evaluation will return None.
+
+        Args:
+            val (Union[ObjParameter, pmo.expression, ObjVariable, np.nan]):
+                Item to be evaluated.
+
+        Returns:
+            Any: If the input is nullable, return np.nan, otherwise evaluate the "value" method of the object.
+        """
+        if pd.isnull(val):
+            return np.nan
+        else:
+            return val.value  # type: ignore
+
+    @staticmethod
+    def _from_pyomo_constraint(
+        val: pmo.constraint, *, eval_body: bool = False
+    ) -> pd.Series:
+        """Evaluate Pyomo constraint object.
+
+        Args:
+            val (pmo.constraint): constraint object to be evaluated
+        Kwargs:
+            eval_body (bool, optional):
+                If True, attempt to evaluate the constraint object `body`, which will evaluate the
+                linear expression contained in the constraint body and produce a numeric value.
+                This will only succeed if the backend model has been successfully optimised,
+                otherwise a string representation of the linear expression will be returned
+                (same as eval_body=False). Defaults to False.
+
+        Returns:
+            pd.Series: Array of upper bound (ub), body, and lower bound (lb).
+        """
+        if pd.isnull(val):
+            vals = [np.nan, np.nan, np.nan]
+        else:
+            if eval_body:
+                try:
+                    body = val.body()
+                except ValueError:
+                    body = val.body.to_string()
+            else:
+                body = val.body.to_string()
+            vals = [val.lb, body, val.ub]
+        return pd.Series(data=vals, index=["lb", "body", "ub"])
+
+    @staticmethod
+    def _from_pyomo_variable_bounds(val: ObjVariable) -> pd.Series:
+        """Evaluate Pyomo decision variable object bounds.
+
+        Args:
+            val (ObjVariable): Variable object to be evaluated.
+
+        Returns:
+            pd.Series: Array of variable upper and lower bound.
+        """
+        if pd.isnull(val):
+            vals = [np.nan, np.nan]
+        else:
+            vals = [val.lb, val.ub]
+        return pd.Series(data=vals, index=["lb", "ub"])
+
+    @staticmethod
+    def _from_pyomo_expr(val: pmo.expression, *, eval_body: bool = False) -> Any:
+        """Evaluate Pyomo expression object.
+
+        Args:
+            val (pmo.expression): expression object to be evaluated
+            eval_body (bool, optional):
+                If True, attempt to evaluate the expression object, which will produce a numeric value.
+                This will only succeed if the backend model has been successfully optimised,
+                otherwise a string representation of the linear expression will be returned
+                (same as eval_body=False). Defaults to False.
+
+        Returns:
+            Any: If the input is nullable, return np.nan, otherwise a numeric value
+            (eval_body=True and problem is optimised) or a string.
+        """
+        if pd.isnull(val):
+            return np.nan
+        else:
+            if eval_body:
+                expr = val()
+                if expr is None:
+                    return val.to_string()
+                else:
+                    return expr
+            else:
+                return val.to_string()
+
+    @contextmanager
+    def _datetime_as_string(self, data: Union[xr.DataArray, xr.Dataset]) -> Iterator:
+        """Context manager to temporarily convert np.dtype("datetime64[ns]") coordinates (e.g. timesteps) to strings with a resolution of minutes.
+
+        Args:
+            data (Union[xr.DataArray, xr.Dataset]): xarray object on whose coordinates the conversion will take place.
+        """
+        datetime_coords = set()
+        for name_, vals_ in data.coords.items():
+            if vals_.dtype.kind == "M":
+                data.coords[name_] = data.coords[name_].dt.strftime("%Y-%m-%d %H:%M")
+                datetime_coords.add(name_)
+        try:
+            yield
+        finally:
+            for name_ in datetime_coords:
+                data.coords[name_] = xr.apply_ufunc(
+                    pd.to_datetime, data.coords[name_], keep_attrs=True
+                )
+
+
+class CoordObj(ABC):
+    """Class with methods to update the `name` property of inheriting classes"""
+
+    def __init__(self) -> None:
+        self._calliope_coords: Optional[Iterable] = None
+
+    def _update_name(self, old_name: str) -> str:
+        """
+        Update string of a list containing a single number with a string of a list containing any arbitrary number of elements
+
+        Args:
+            old_name (str): String representation of a list containing a single number
+
+        Returns:
+            str:
+                If `self.calliope_coords` is None, returns `old_name`.
+                Otherwise returns string representation of a list containing the contents of `self.calliope_coords`
+        """
+        if self._calliope_coords is None:
+            return old_name
+
+        if not self._calliope_coords:  # empty list = dimensionless component
+            coord_list = ""
+        else:
+            coord_list = f"[{', '.join(str(i) for i in self._calliope_coords)}]"
+        return re.sub(r"\[\d+\]", coord_list, old_name)
+
+    @property
+    def calliope_coords(self):
+        return self._calliope_coords
+
+    @calliope_coords.setter
+    def calliope_coords(self, val):
+        self._calliope_coords = val
+
+
+class ObjParameter(pmo.parameter, CoordObj):
+    """
+    A pyomo parameter (`a object for storing a mutable, numeric value that can be used to build a symbolic expression`)
+    with added `dtype` property and a `name` property setter (via the `pmo.parameter.getname` method) which replaces a list position as a name with a list of strings.
+    """
+
+    def __init__(self, value, **kwds):
+        assert not pd.isnull(value)
+        pmo.parameter.__init__(self, value, **kwds)
+        CoordObj.__init__(self)
+
+    @property
+    def dtype(self):
+        return "O"
+
+    def getname(self, *args, **kwargs):
+        return self._update_name(pmo.parameter.getname(self, *args, **kwargs))
+
+
+class ObjVariable(pmo.variable, CoordObj):
+    """
+    A pyomo variable with a `name` property setter (via the `pmo.variable.getname` method) which replaces a list position as a name with a list of strings.
+
+    """
+
+    def __init__(self, **kwds):
+        pmo.variable.__init__(self, **kwds)
+        CoordObj.__init__(self)
+
+    def getname(self, *args, **kwargs):
+        return self._update_name(pmo.variable.getname(self, *args, **kwargs))
diff --git a/src/calliope/cli.py b/src/calliope/cli.py
index 775f614d2..dcf4eeb93 100644
--- a/src/calliope/cli.py
+++ b/src/calliope/cli.py
@@ -284,7 +284,7 @@ def run(
                     bold=True,
                 )
             model.build()
-            model.to_lp(save_lp)
+            model.backend.to_lp(save_lp)
             print_end_time(start_time)
 
         # Else run the model, then save outputs
diff --git a/src/calliope/core/model.py b/src/calliope/core/model.py
index 49a40be48..97e976129 100644
--- a/src/calliope/core/model.py
+++ b/src/calliope/core/model.py
@@ -19,11 +19,13 @@
 
 import calliope
 from calliope import exceptions
-from calliope.backend import backends, latex_backend, parsing
+from calliope.backend import parsing
+from calliope.backend.latex_backend_model import LatexBackendModel, MathDocumentation
+from calliope.backend.pyomo_backend_model import PyomoBackendModel
 from calliope.core import io
 from calliope.core.attrdict import AttrDict
 from calliope.core.util.logging import log_time
-from calliope.core.util.tools import copy_docstring, relative_path, validate_dict
+from calliope.core.util.tools import relative_path, validate_dict
 from calliope.postprocess import results as postprocess_results
 from calliope.preprocess import model_run_from_dict, model_run_from_yaml
 from calliope.preprocess.model_data import ModelDataFactory
@@ -31,7 +33,8 @@
 logger = logging.getLogger(__name__)
 
 T = TypeVar(
-    "T", bound=Union[backends.PyomoBackendModel, latex_backend.LatexBackendModel]
+    "T",
+    bound=Union[PyomoBackendModel, LatexBackendModel],
 )
 
 
@@ -48,7 +51,7 @@ class Model(object):
     A Calliope Model.
     """
 
-    _BACKENDS: dict[str, Callable] = {"pyomo": backends.PyomoBackendModel}
+    _BACKENDS: dict[str, Callable] = {"pyomo": PyomoBackendModel}
     _MATH_SCHEMA = AttrDict.from_yaml(
         Path(calliope.__file__).parent / "config" / "math_schema.yaml"
     )
@@ -96,7 +99,7 @@ def __init__(
         self.run_config: AttrDict
         self.math: AttrDict
         self._config_path: Optional[str]
-        self.math_documentation = latex_backend.MathDocumentation(self._build)
+        self.math_documentation = MathDocumentation(self._build)
 
         # try to set logging output format assuming python interactive. Will
         # use CLI logging format if model called from CLI
@@ -171,6 +174,8 @@ def _init_from_model_run(
         self._add_observed_dict("math", math)
 
         self.inputs = self._model_data.filter_by_attrs(is_result=0)
+        self.math_documentation.inputs = self._model_data
+
         log_time(
             logger,
             self._timings,
@@ -221,6 +226,8 @@ def _add_model_data_methods(self):
 
         self.inputs = self._model_data.filter_by_attrs(is_result=0)
         results = self._model_data.filter_by_attrs(is_result=1)
+        self.math_documentation.inputs = self._model_data
+
         if len(results.data_vars) > 0:
             self.results = results
         log_time(
@@ -363,11 +370,10 @@ def build(
                 "This model object already has a built optimisation problem. Use model.build(force=True) "
                 "to force the existing optimisation problem to be overwritten with a new one."
             )
-        backend = self._BACKENDS[backend_interface]()
+        backend = self._BACKENDS[backend_interface](self._model_data)
         self.backend = self._build(backend)
 
     def _build(self, backend: T) -> T:
-        backend.add_all_parameters(self._model_data, self.run_config)
         log_time(
             logger,
             self._timings,
@@ -387,13 +393,9 @@ def _build(self, backend: T) -> T:
             component = components.removesuffix("s")
             if components in ["variables", "global_expressions"]:
                 backend.valid_math_element_names.update(self.math[components].keys())
-            for name, dict_ in self.math[components].items():
-                if dict_.get("active", True):
-                    getattr(backend, f"add_{component}")(self._model_data, name, dict_)
-                else:
-                    logger.debug(
-                        f"({component}, {name}): Component deactivated and therefore not built."
-                    )
+            for name in self.math[components]:
+                getattr(backend, f"add_{component}")(name)
+
             log_time(
                 logger,
                 self._timings,
@@ -402,14 +404,6 @@ def _build(self, backend: T) -> T:
             )
         return backend
 
-    @copy_docstring(backends.BackendModel.verbose_strings)
-    def verbose_strings(self) -> None:
-        if not hasattr(self, "backend"):
-            raise NotImplementedError(
-                "Call `build()` to generate an optimisation problem before calling this function."
-            )
-        self.backend.verbose_strings()
-
     def solve(self, force: bool = False, warmstart: bool = False) -> None:
         """
         Run the built optimisation problem.
@@ -559,19 +553,6 @@ def to_csv(self, path, dropna=True):
         """
         io.save_csv(self._model_data, path, dropna)
 
-    @copy_docstring(backends.BackendModel.to_lp)
-    def to_lp(self, path: Union[str, Path]) -> None:
-        """
-        Raises:
-            exceptions.ModelError: This method cannot be called prior to calling `build()`.
-        """
-
-        if not hasattr(self, "backend"):
-            raise exceptions.ModelError(
-                "Build the optimisation problem by calling `build()` before trying to generate an LP file."
-            )
-        self.backend.to_lp(path)
-
     def info(self) -> str:
         """Generate basic description of the model, combining its name and a rough indication of the model size.
 
diff --git a/src/calliope/core/util/logging.py b/src/calliope/core/util/logging.py
index 8a4e134aa..7d724ac6e 100644
--- a/src/calliope/core/util/logging.py
+++ b/src/calliope/core/util/logging.py
@@ -62,7 +62,7 @@ def set_log_verbosity(
         interactively.
 
     """
-    backend_logger = logging.getLogger("calliope.backend.backends")
+    backend_logger = logging.getLogger("calliope.backend.backend_model")
     if include_solver_output is True:
         backend_logger.setLevel("DEBUG")
     else:
diff --git a/src/calliope/core/util/tools.py b/src/calliope/core/util/tools.py
index c8f861536..ef4223a99 100644
--- a/src/calliope/core/util/tools.py
+++ b/src/calliope/core/util/tools.py
@@ -4,10 +4,9 @@
 
 import functools
 import importlib
-import operator
 import os
 import sys
-from typing import Callable, TypeVar
+from typing import TypeVar
 
 import jsonschema
 from typing_extensions import ParamSpec
@@ -18,14 +17,6 @@
 T = TypeVar("T")
 
 
-def get_from_dict(data_dict, map_list):
-    return functools.reduce(operator.getitem, map_list, data_dict)
-
-
-def apply_to_dict(data_dict, map_list, func, args):
-    getattr(get_from_dict(data_dict, map_list[:-1])[map_list[-1]], func)(*args)
-
-
 memoize = functools.lru_cache(maxsize=2048)
 
 
@@ -111,26 +102,6 @@ def plugin_load(name, builtin_module):
     return func
 
 
-def copy_docstring(wrapper: Callable[P, T]):
-    """
-    Decorator to copy across a function docstring to the wrapped function.
-    Any additional documentation in the wrapped function will be appended to the copied
-    docstring.
-    """
-
-    def decorator(func: Callable) -> Callable[P, T]:
-        func_doc = ""
-        if wrapper.__doc__ is not None:
-            func_doc += wrapper.__doc__
-        if func.__doc__ is not None:
-            func_doc += func.__doc__
-        func.__doc__ = func_doc
-
-        return func
-
-    return decorator
-
-
 def validate_dict(to_validate: dict, schema: dict, dict_descriptor: str) -> None:
     """
     Validate a dictionary under a given schema.


Test patch:
diff --git a/tests/common/lp_files/balance_conversion.lp b/tests/common/lp_files/balance_conversion.lp
index 55c804eab..c179e0bd3 100644
--- a/tests/common/lp_files/balance_conversion.lp
+++ b/tests/common/lp_files/balance_conversion.lp
@@ -1,7 +1,7 @@
 \* Source Pyomo model name=None *\
 
 min
-objectives(0):
+objectives(dummy_obj)(0):
 +2.0 ONE_VAR_CONSTANT
 
 s.t.
diff --git a/tests/common/lp_files/carrier_production_max.lp b/tests/common/lp_files/carrier_production_max.lp
index b04ae1aa9..f885f58db 100644
--- a/tests/common/lp_files/carrier_production_max.lp
+++ b/tests/common/lp_files/carrier_production_max.lp
@@ -1,7 +1,7 @@
 \* Source Pyomo model name=None *\
 
 min
-objectives(0):
+objectives(dummy_obj)(0):
 +2.0 ONE_VAR_CONSTANT
 
 s.t.
diff --git a/tests/common/lp_files/energy_cap.lp b/tests/common/lp_files/energy_cap.lp
index 4654af5ad..280fbb0ee 100644
--- a/tests/common/lp_files/energy_cap.lp
+++ b/tests/common/lp_files/energy_cap.lp
@@ -1,7 +1,7 @@
 \* Source Pyomo model name=None *\
 
 min
-objectives(0):
+objectives(foo)(0):
 +1 variables(energy_cap)(a__test_supply_elec)
 +1 variables(energy_cap)(b__test_supply_elec)
 
diff --git a/tests/common/lp_files/resource_max.lp b/tests/common/lp_files/resource_max.lp
index 0bda39057..96b633772 100644
--- a/tests/common/lp_files/resource_max.lp
+++ b/tests/common/lp_files/resource_max.lp
@@ -1,7 +1,7 @@
 \* Source Pyomo model name=None *\
 
 min
-objectives(0):
+objectives(dummy_obj)(0):
 +2.0 ONE_VAR_CONSTANT
 
 s.t.
diff --git a/tests/common/lp_files/storage_max.lp b/tests/common/lp_files/storage_max.lp
index b23f6c6cc..8a078cc54 100644
--- a/tests/common/lp_files/storage_max.lp
+++ b/tests/common/lp_files/storage_max.lp
@@ -1,7 +1,7 @@
 \* Source Pyomo model name=None *\
 
 min
-objectives(0):
+objectives(dummy_obj)(0):
 +2 ONE_VAR_CONSTANT
 
 s.t.
diff --git a/tests/common/util.py b/tests/common/util.py
index eefce6e5c..51cc4bfdc 100644
--- a/tests/common/util.py
+++ b/tests/common/util.py
@@ -102,26 +102,27 @@ def build_lp(
         math (Optional[dict], optional): All constraint/global expression/objective math to apply. Defaults to None.
         backend (Literal["pyomo"], optional): Backend to use to create the LP file. Defaults to "pyomo".
     """
-    backend_instance = model._BACKENDS[backend]()
-    backend_instance.add_all_parameters(model.inputs, model.run_config)
+    backend_instance = model._BACKENDS[backend](model._model_data)
     for name, dict_ in model.math["variables"].items():
-        backend_instance.add_variable(model.inputs, name, dict_)
+        backend_instance.add_variable(name, dict_)
 
     if math is not None:
         for component_group, component_math in math.items():
             for name, dict_ in component_math.items():
                 getattr(backend_instance, f"add_{component_group.removesuffix('s')}")(
-                    model.inputs, name, dict_
+                    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(
-            model.inputs,
             "dummy_obj",
             {"equations": [{"expression": "1 + 1"}], "sense": "minimize"},
         )
-    backend_instance._instance.objectives[0].activate()
+        backend_instance._instance.objectives["dummy_obj"][0].activate()
+    elif "objectives" in math.keys():
+        objective = list(math["objectives"].keys())[0]
+        backend_instance._instance.objectives[objective][0].activate()
 
     backend_instance.verbose_strings()
 
diff --git a/tests/conftest.py b/tests/conftest.py
index aa9092abd..326fc3aa3 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -4,7 +4,7 @@
 import pytest
 import xarray as xr
 from calliope import AttrDict
-from calliope.backend import backends, latex_backend
+from calliope.backend import latex_backend_model, pyomo_backend_model
 
 from .common.util import build_test_model as build_model
 
@@ -228,10 +228,8 @@ def dummy_model_data():
     return model_data
 
 
-def populate_backend_model(backend, model_data):
-    backend.add_all_parameters(model_data, model_data.run_config)
+def populate_backend_model(backend):
     backend.add_variable(
-        model_data,
         "multi_dim_var",
         {
             "foreach": ["nodes", "techs"],
@@ -239,19 +237,17 @@ def populate_backend_model(backend, model_data):
             "bounds": {"min": "0", "max": ".inf"},
         },
     )
-    backend.add_variable(
-        model_data, "no_dim_var", {"bounds": {"min": "-1", "max": "1"}}
-    )
+    backend.add_variable("no_dim_var", {"bounds": {"min": "-1", "max": "1"}})
     return backend
 
 
 @pytest.fixture(scope="module")
 def dummy_pyomo_backend_model(dummy_model_data):
-    backend = backends.PyomoBackendModel()
-    return populate_backend_model(backend, dummy_model_data)
+    backend = pyomo_backend_model.PyomoBackendModel(dummy_model_data)
+    return populate_backend_model(backend)
 
 
 @pytest.fixture(scope="module")
 def dummy_latex_backend_model(dummy_model_data):
-    backend = latex_backend.LatexBackendModel()
-    return populate_backend_model(backend, dummy_model_data)
+    backend = latex_backend_model.LatexBackendModel(dummy_model_data)
+    return populate_backend_model(backend)
diff --git a/tests/test_backend_latex_backend.py b/tests/test_backend_latex_backend.py
index 0f87cd6dd..7de4779ac 100644
--- a/tests/test_backend_latex_backend.py
+++ b/tests/test_backend_latex_backend.py
@@ -4,7 +4,7 @@
 import pytest
 import xarray as xr
 from calliope import exceptions
-from calliope.backend import latex_backend
+from calliope.backend import latex_backend_model
 
 from .common.util import build_test_model, check_error_or_warning
 
@@ -67,8 +67,8 @@ def test_invalid_format(self, build_all, tmpdir_factory, filepath, format):
 
 class TestLatexBackendModel:
     @pytest.fixture(scope="class")
-    def valid_latex_backend(self):
-        return latex_backend.LatexBackendModel(include="valid")
+    def valid_latex_backend(self, dummy_model_data):
+        return latex_backend_model.LatexBackendModel(dummy_model_data, include="valid")
 
     @pytest.mark.parametrize(
         "backend_obj", ["valid_latex_backend", "dummy_latex_backend_model"]
@@ -85,7 +85,6 @@ def test_add_parameter(self, request, backend_obj):
     def test_add_variable(self, request, dummy_model_data, backend_obj):
         latex_backend_model = request.getfixturevalue(backend_obj)
         latex_backend_model.add_variable(
-            dummy_model_data,
             "var",
             {
                 "foreach": ["nodes", "techs"],
@@ -101,9 +100,8 @@ def test_add_variable(self, request, dummy_model_data, backend_obj):
         assert "var" in latex_backend_model.valid_math_element_names
         assert "math_string" in latex_backend_model.variables["var"].attrs
 
-    def test_add_variable_not_valid(self, dummy_model_data, valid_latex_backend):
+    def test_add_variable_not_valid(self, valid_latex_backend):
         valid_latex_backend.add_variable(
-            dummy_model_data,
             "invalid_var",
             {
                 "foreach": ["nodes", "techs"],
@@ -122,7 +120,6 @@ def test_add_variable_not_valid(self, dummy_model_data, valid_latex_backend):
     def test_add_expression(self, request, dummy_model_data, backend_obj):
         latex_backend_model = request.getfixturevalue(backend_obj)
         latex_backend_model.add_global_expression(
-            dummy_model_data,
             "expr",
             {
                 "foreach": ["nodes", "techs"],
@@ -144,7 +141,6 @@ def test_add_expression(self, request, dummy_model_data, backend_obj):
     def test_add_constraint(self, request, dummy_model_data, backend_obj):
         latex_backend_model = request.getfixturevalue(backend_obj)
         latex_backend_model.add_constraint(
-            dummy_model_data,
             "constr",
             {
                 "foreach": ["nodes", "techs"],
@@ -160,9 +156,8 @@ def test_add_constraint(self, request, dummy_model_data, backend_obj):
         assert "constr" not in latex_backend_model.valid_math_element_names
         assert "math_string" in latex_backend_model.constraints["constr"].attrs
 
-    def test_add_constraint_not_valid(self, dummy_model_data, valid_latex_backend):
+    def test_add_constraint_not_valid(self, valid_latex_backend):
         valid_latex_backend.add_constraint(
-            dummy_model_data,
             "invalid_constr",
             {
                 "foreach": ["nodes", "techs"],
@@ -173,14 +168,13 @@ def test_add_constraint_not_valid(self, dummy_model_data, valid_latex_backend):
                 ],
             },
         )
-        assert not valid_latex_backend.constraints["invalid_constr"].any()
+        assert valid_latex_backend.constraints["invalid_constr"].isnull().all()
         assert (
             "math_string" not in valid_latex_backend.constraints["invalid_constr"].attrs
         )
 
-    def test_add_constraint_one_not_valid(self, dummy_model_data, valid_latex_backend):
+    def test_add_constraint_one_not_valid(self, valid_latex_backend):
         valid_latex_backend.add_constraint(
-            dummy_model_data,
             "valid_constr",
             {
                 "foreach": ["nodes", "techs"],
@@ -196,9 +190,8 @@ def test_add_constraint_one_not_valid(self, dummy_model_data, valid_latex_backen
             not in valid_latex_backend.constraints["valid_constr"].attrs["math_string"]
         )
 
-    def test_add_objective(self, dummy_model_data, dummy_latex_backend_model):
+    def test_add_objective(self, dummy_latex_backend_model):
         dummy_latex_backend_model.add_objective(
-            dummy_model_data,
             "obj",
             {
                 "equations": [{"expression": "sum(var, over=[nodes, techs])"}],
@@ -209,29 +202,8 @@ def test_add_objective(self, dummy_model_data, dummy_latex_backend_model):
         assert "obj" not in dummy_latex_backend_model.valid_math_element_names
         assert len(dummy_latex_backend_model.objectives.data_vars) == 1
 
-    def test_get_parameter(self, dummy_latex_backend_model):
-        param = dummy_latex_backend_model.get_parameter("param")
-        assert param.equals(xr.DataArray(1))
-
     def test_create_obj_list(self, dummy_latex_backend_model):
-        assert dummy_latex_backend_model.create_obj_list("var", "variables") is None
-
-    def test_get_constraint(self, dummy_latex_backend_model):
-        constr = dummy_latex_backend_model.get_constraint("constr")
-        assert constr.equals(dummy_latex_backend_model.constraints["constr"])
-
-    def test_get_variable(self, dummy_latex_backend_model):
-        var = dummy_latex_backend_model.get_variable("var")
-        assert var.equals(dummy_latex_backend_model.variables["var"])
-
-    def test_get_expression(self, dummy_latex_backend_model):
-        expr = dummy_latex_backend_model.get_global_expression("expr")
-        assert expr.equals(dummy_latex_backend_model.global_expressions["expr"])
-
-    def test_solve(self, dummy_latex_backend_model):
-        with pytest.raises(exceptions.BackendError) as excinfo:
-            dummy_latex_backend_model.solve("cbc")
-        assert check_error_or_warning(excinfo, "Cannot solve a LaTex backend model")
+        assert dummy_latex_backend_model._create_obj_list("var", "variables") is None
 
     @pytest.mark.parametrize(
         ["format", "expected"],
@@ -321,13 +293,11 @@ def test_solve(self, dummy_latex_backend_model):
         ],
     )
     def test_generate_math_doc(self, dummy_model_data, format, expected):
-        latex_backend_model = latex_backend.LatexBackendModel()
-        latex_backend_model.add_global_expression(
-            dummy_model_data,
-            "expr",
-            {"equations": [{"expression": "1 + 2"}], "description": "foobar"},
+        backend_model = latex_backend_model.LatexBackendModel(dummy_model_data)
+        backend_model.add_global_expression(
+            "expr", {"equations": [{"expression": "1 + 2"}], "description": "foobar"}
         )
-        doc = latex_backend_model.generate_math_doc(format=format)
+        doc = backend_model.generate_math_doc(format=format)
         assert doc == expected
 
     @pytest.mark.parametrize(
@@ -403,7 +373,7 @@ def test_generate_math_doc(self, dummy_model_data, format, expected):
     )
     def test_generate_math_string(self, dummy_latex_backend_model, kwargs, expected):
         da = xr.DataArray()
-        dummy_latex_backend_model._generate_math_string(da, **kwargs)
+        dummy_latex_backend_model._generate_math_string(None, da, **kwargs)
         assert da.math_string == expected
 
     @pytest.mark.parametrize(
@@ -418,11 +388,9 @@ def test_render(self, dummy_latex_backend_model, instring, kwargs, expected):
         rendered = dummy_latex_backend_model._render(instring, **kwargs)
         assert rendered == expected
 
-    def test_get_capacity_bounds(self, dummy_latex_backend_model, dummy_model_data):
+    def test_get_capacity_bounds(self, dummy_latex_backend_model):
         bounds = {"min": 1, "max": 2e6}
-        lb, ub = dummy_latex_backend_model._get_capacity_bounds(
-            bounds, "multi_dim_var", dummy_model_data
-        )
+        lb, ub = dummy_latex_backend_model._get_capacity_bounds("multi_dim_var", bounds)
         assert lb == {"expression": r"1 \leq \textbf{multi_dim_var}_\text{node,tech}"}
         assert ub == {
             "expression": r"\textbf{multi_dim_var}_\text{node,tech} \leq 2\mathord{\times}10^{+06}"
diff --git a/tests/test_backend_parsing.py b/tests/test_backend_parsing.py
index 7f87040a5..acd351183 100644
--- a/tests/test_backend_parsing.py
+++ b/tests/test_backend_parsing.py
@@ -5,7 +5,7 @@
 import pyparsing as pp
 import pytest
 import ruamel.yaml as yaml
-from calliope.backend import backends, expression_parser, parsing, where_parser
+from calliope.backend import backend_model, expression_parser, parsing, where_parser
 
 from .common.util import check_error_or_warning
 
@@ -212,12 +212,14 @@ def _equation_slice_obj(name):
 
 @pytest.fixture
 def dummy_backend_interface(dummy_model_data):
-    # ignore the need to define the abstract methods from backends.BackendModel
-    with patch.multiple(backends.BackendModel, __abstractmethods__=set()):
+    # ignore the need to define the abstract methods from backend_model.BackendModel
+    with patch.multiple(backend_model.BackendModel, __abstractmethods__=set()):
 
-        class DummyBackendModel(backends.BackendModel):
+        class DummyBackendModel(backend_model.BackendModel):
             def __init__(self):
-                backends.BackendModel.__init__(self, instance=None)
+                backend_model.BackendModel.__init__(
+                    self, dummy_model_data, instance=None
+                )
 
                 self._dataset = dummy_model_data.copy(deep=True)
                 self._dataset["with_inf"] = self._dataset["with_inf"].fillna(
diff --git a/tests/test_backend_pyomo.py b/tests/test_backend_pyomo.py
index 84707afc2..094841826 100755
--- a/tests/test_backend_pyomo.py
+++ b/tests/test_backend_pyomo.py
@@ -3,6 +3,7 @@
 
 import calliope.exceptions as exceptions
 import numpy as np
+import pandas as pd
 import pyomo.kernel as pmo
 import pytest  # noqa: F401
 import xarray as xr
@@ -1802,6 +1803,8 @@ def test_no_duplicate_log_message(self, caplog, gurobi_model):
 
 
 class TestNewBackend:
+    LOGGER = logging.getLogger("calliope.backend.backend_model")
+
     @pytest.fixture(scope="class")
     def simple_supply_longnames(self):
         m = build_model({}, "simple_supply,two_hours,investment_costs")
@@ -1820,8 +1823,8 @@ def test_new_build_optimal(self, simple_supply):
         "component_type", ["variable", "global_expression", "parameter", "constraint"]
     )
     def test_new_build_get_missing_component(self, simple_supply, component_type):
-        returned_ = getattr(simple_supply.backend, f"get_{component_type}")("foo")
-        assert returned_ is None
+        with pytest.raises(KeyError):
+            getattr(simple_supply.backend, f"get_{component_type}")("foo")
 
     def test_new_build_get_variable(self, simple_supply):
         var = simple_supply.backend.get_variable("energy_cap")
@@ -1981,14 +1984,12 @@ def test_solve_warmstart_not_possible(self, simple_supply):
         assert check_error_or_warning(excinfo, "cbc, does not support warmstart")
 
     def test_solve_non_optimal(self, simple_supply):
-        def _update_param(param):
-            param.value = param.value * 1000
-
-        simple_supply.backend.apply_func(
-            _update_param,
-            simple_supply.backend.parameters.resource.loc[
-                {"techs": "test_demand_elec"}
-            ],
+        simple_supply.backend.update_parameter(
+            "resource",
+            simple_supply.inputs.resource.where(
+                simple_supply.inputs.techs == "test_demand_elec"
+            )
+            * 1000,
         )
         with pytest.warns(exceptions.BackendWarning) as excinfo:
             simple_supply.solve(force=True)
@@ -2033,7 +2034,6 @@ def test_raise_error_on_constraint_with_nan(self, simple_supply):
         constraint_name = "constraint-without-nan"
 
         simple_supply.backend.add_constraint(
-            simple_supply.inputs,
             constraint_name,
             constraint_dict,
         )
@@ -2055,7 +2055,6 @@ def test_raise_error_on_constraint_with_nan(self, simple_supply):
 
         with pytest.raises(exceptions.BackendError) as error:
             simple_supply.backend.add_constraint(
-                simple_supply.inputs,
                 constraint_name,
                 constraint_dict,
             )
@@ -2081,7 +2080,6 @@ def test_raise_error_on_expression_with_nan(self, simple_supply):
 
         # add expression with nan
         simple_supply.backend.add_global_expression(
-            simple_supply.inputs,
             expression_name,
             expression_dict,
         )
@@ -2100,7 +2098,6 @@ def test_raise_error_on_expression_with_nan(self, simple_supply):
 
         with pytest.raises(exceptions.BackendError) as error:
             simple_supply.backend.add_global_expression(
-                simple_supply.inputs,
                 expression_name,
                 expression_dict,
             )
@@ -2127,7 +2124,6 @@ def test_raise_error_on_excess_dimensions(self, simple_supply):
         constraint_name = "constraint-without-excess-dimensions"
 
         simple_supply.backend.add_constraint(
-            simple_supply.inputs,
             constraint_name,
             constraint_dict,
         )
@@ -2148,7 +2144,6 @@ def test_raise_error_on_excess_dimensions(self, simple_supply):
 
         with pytest.raises(exceptions.BackendError) as error:
             simple_supply.backend.add_constraint(
-                simple_supply.inputs,
                 constraint_name,
                 constraint_dict,
             )
@@ -2163,11 +2158,12 @@ def test_raise_error_on_excess_dimensions(self, simple_supply):
     )
     def test_create_and_delete_pyomo_list(self, simple_supply, component):
         backend_instance = simple_supply.backend._instance
-        simple_supply.backend.create_obj_list("foo", component)
+        simple_supply.backend._create_obj_list("foo", component)
         assert "foo" in getattr(backend_instance, component).keys()
 
-        simple_supply.backend.delete_obj_list("foo", component)
+        simple_supply.backend.delete_component("foo", component)
         assert "foo" not in getattr(backend_instance, component).keys()
+        assert "foo" not in getattr(simple_supply.backend, component).keys()
 
     @pytest.mark.parametrize(
         "component", ["parameters", "variables", "global_expressions", "constraints"]
@@ -2175,7 +2171,7 @@ def test_create_and_delete_pyomo_list(self, simple_supply, component):
     def test_delete_inexistent_pyomo_list(self, simple_supply, component):
         backend_instance = simple_supply.backend._instance
         assert "bar" not in getattr(backend_instance, component).keys()
-        simple_supply.backend.delete_obj_list("bar", component)
+        simple_supply.backend.delete_component("bar", component)
         assert "bar" not in getattr(backend_instance, component).keys()
 
     @pytest.mark.parametrize(
@@ -2189,7 +2185,7 @@ def test_add_allnull_expr_or_constr(self, simple_supply, component, eq):
             "where": "True",
             "equations": [{"expression": eq, "where": "False"}],
         }
-        adder(simple_supply._model_data, "foo", constr_dict)
+        adder("foo", constr_dict)
 
         assert "foo" not in getattr(simple_supply.backend._instance, component).keys()
         assert "foo" not in simple_supply.backend._dataset.data_vars.keys()
@@ -2213,7 +2209,6 @@ def test_add_allnull_param_with_shape(self, simple_supply):
 
     def test_add_allnull_var(self, simple_supply):
         simple_supply.backend.add_variable(
-            simple_supply._model_data,
             "foo",
             {"foreach": ["nodes"], "where": "False"},
         )
@@ -2223,9 +2218,7 @@ def test_add_allnull_var(self, simple_supply):
     def test_add_allnull_obj(self, simple_supply):
         eq = {"expression": "bigM", "where": "False"}
         simple_supply.backend.add_objective(
-            simple_supply._model_data,
-            "foo",
-            {"equations": [eq, eq], "sense": "minimise"},
+            "foo", {"equations": [eq, eq], "sense": "minimise"}
         )
         assert len(simple_supply.backend._instance.objectives) == 1
         assert "foo" not in simple_supply.backend._dataset.data_vars.keys()
@@ -2234,21 +2227,17 @@ def test_add_two_same_obj(self, simple_supply):
         eq = {"expression": "bigM", "where": "True"}
         with pytest.raises(exceptions.BackendError) as excinfo:
             simple_supply.backend.add_objective(
-                simple_supply._model_data,
-                "foo",
-                {"equations": [eq, eq], "sense": "minimise"},
+                "foo", {"equations": [eq, eq], "sense": "minimise"}
             )
         assert check_error_or_warning(
             excinfo,
-            "More than one foo objective is valid for this optimisation problem; only one is allowed.",
+            "(objective, foo:1): trying to set two equations for the same component.",
         )
 
     def test_add_valid_obj(self, simple_supply):
         eq = {"expression": "bigM", "where": "True"}
         simple_supply.backend.add_objective(
-            simple_supply._model_data,
-            "foo",
-            {"equations": [eq], "sense": "minimise"},
+            "foo", {"equations": [eq], "sense": "minimise"}
         )
         assert "foo" in simple_supply.backend.objectives
         assert not simple_supply.backend.objectives.foo.item().active
@@ -2332,3 +2321,182 @@ def test_verbose_strings_no_len(self, simple_supply_longnames):
 
         assert obj.item().name == "parameters[annualisation_weight]"
         assert obj.coords_in_name
+
+    def test_update_parameter(self, simple_supply):
+        updated_param = simple_supply.inputs.energy_eff * 1000
+        simple_supply.backend.update_parameter("energy_eff", updated_param)
+
+        expected = simple_supply.backend.get_parameter(
+            "energy_eff", as_backend_objs=False
+        )
+        assert expected.where(updated_param.notnull()).equals(updated_param)
+
+    def test_update_parameter_one_val(self, caplog, simple_supply):
+        updated_param = 1000
+        new_dims = {"nodes", "techs"}
+        caplog.set_level(logging.WARNING)
+
+        simple_supply.backend.update_parameter("energy_eff", updated_param)
+
+        assert (
+            f"New values will be broadcast along the {new_dims} dimension(s)"
+            in caplog.text
+        )
+        expected = simple_supply.backend.get_parameter(
+            "energy_eff", as_backend_objs=False
+        )
+        assert (expected == updated_param).all()
+
+    def test_update_parameter_replace_defaults(self, simple_supply):
+        updated_param = simple_supply.inputs.energy_eff.fillna(0.1)
+
+        simple_supply.backend.update_parameter("energy_eff", updated_param)
+
+        expected = simple_supply.backend.get_parameter(
+            "energy_eff", as_backend_objs=False
+        )
+        assert expected.equals(updated_param)
+
+    def test_update_parameter_add_dim(self, caplog, simple_supply):
+        """energy_eff doesn't have the time dimension in the simple model, we add it here."""
+        updated_param = simple_supply.inputs.energy_eff.where(
+            simple_supply.inputs.timesteps.notnull()
+        )
+
+        refs_to_update = {"balance_demand", "balance_transmission"}
+        caplog.set_level(logging.WARNING)
+
+        simple_supply.backend.update_parameter("energy_eff", updated_param)
+
+        assert (
+            f"Defining values for a previously fully/partially undefined parameter. The optimisation problem components {refs_to_update} will be re-built."
+            in caplog.text
+        )
+
+        expected = simple_supply.backend.get_parameter(
+            "energy_eff", as_backend_objs=False
+        )
+        assert "timesteps" in expected.dims
+
+    def test_update_parameter_replace_undefined(self, caplog, simple_supply):
+        """parasitic_eff isn't defined in the inputs, so is a dimensionless value in the pyomo object, assigned its default value"""
+        updated_param = simple_supply.inputs.energy_eff
+
+        refs_to_update = {"carrier_production_max"}
+        caplog.set_level(logging.WARNING)
+
+        simple_supply.backend.update_parameter("parasitic_eff", updated_param)
+
+        assert (
+            f"Defining values for a previously fully/partially undefined parameter. The optimisation problem components {refs_to_update} will be re-built."
+            in caplog.text
+        )
+
+        expected = simple_supply.backend.get_parameter(
+            "parasitic_eff", as_backend_objs=False
+        )
+        default_val = simple_supply._model_data.defaults["parasitic_eff"]
+        assert expected.equals(updated_param.fillna(default_val))
+
+    def test_update_parameter_no_refs_to_update(self, simple_supply):
+        """units_equals isn't defined in the inputs, so is a dimensionless value in the pyomo object, assigned its default value.
+
+        Updating it doesn't change the model in any way, because none of the existing constraints/expressions depend on it. Therefore, no warning is raised
+        """
+
+        updated_param = 1
+
+        simple_supply.backend.update_parameter("units_equals", updated_param)
+
+        expected = simple_supply.backend.get_parameter(
+            "units_equals", as_backend_objs=False
+        )
+        assert expected == 1
+
+    @pytest.mark.parametrize("bound", ["min", "max"])
+    def test_update_variable_single_bound_single_val(self, simple_supply, bound):
+        translator = {"min": "lb", "max": "ub"}
+
+        simple_supply.backend.update_variable_bounds("carrier_prod", **{bound: 1})
+
+        bound_vals = simple_supply.backend.get_variable_bounds("carrier_prod")[
+            translator[bound]
+        ]
+
+        assert (bound_vals == 1).where(bound_vals.notnull()).all()
+
+    def test_update_variable_bounds_single_val(self, simple_supply):
+        simple_supply.backend.update_variable_bounds("carrier_prod", min=2, max=2)
+        bound_vals = simple_supply.backend.get_variable_bounds("carrier_prod")
+        assert (bound_vals == 2).where(bound_vals.notnull()).all().all()
+
+    def test_update_variable_single_bound_multi_val(self, caplog, simple_supply):
+        bound_array = simple_supply.inputs.resource.sel(techs="test_demand_elec")
+        simple_supply.backend.update_variable_bounds("carrier_con", min=bound_array)
+        bound_vals = simple_supply.backend.get_variable_bounds("carrier_con").lb
+        assert (
+            "New `min` bounds for variable `carrier_con` will be broadcast"
+            in caplog.text
+        )
+        assert bound_vals.equals(
+            bound_array.where(bound_vals.notnull()).transpose(*bound_vals.dims)
+        )
+
+    def test_update_variable_error_update_parameter_instead(self, simple_supply):
+        with pytest.raises(exceptions.BackendError) as excinfo:
+            simple_supply.backend.update_variable_bounds("energy_cap", min=1)
+        assert check_error_or_warning(
+            excinfo,
+            "Cannot update variable bounds that have been set by parameters. use `update_parameter('energy_cap_min')` to update the min bound of energy_cap",
+        )
+
+    @staticmethod
+    def _is_fixed(val):
+        if pd.notnull(val):
+            return val.fixed
+
+    def test_fix_variable(self, simple_supply):
+        simple_supply.backend.fix_variable("energy_cap")
+        fixed = simple_supply.backend._apply_func(
+            self._is_fixed, simple_supply.backend.variables.energy_cap
+        )
+        simple_supply.backend.unfix_variable("energy_cap")  # reset
+        assert fixed.where(fixed.notnull()).all()
+
+    def test_fix_variable_where(self, simple_supply):
+        where = simple_supply.inputs.energy_cap_max.notnull()
+        simple_supply.backend.fix_variable("energy_cap", where=where)
+        fixed = simple_supply.backend._apply_func(
+            self._is_fixed, simple_supply.backend.variables.energy_cap
+        )
+        simple_supply.backend.unfix_variable("energy_cap")  # reset
+        assert not fixed.sel(techs="test_demand_elec").any()
+        assert fixed.where(where).all()
+
+    def test_fix_variable_before_solve(self, simple_supply_longnames):
+        with pytest.raises(exceptions.BackendError) as excinfo:
+            simple_supply_longnames.backend.fix_variable("energy_cap")
+
+        assert check_error_or_warning(
+            excinfo,
+            "Cannot fix variable values without already having solved the model successfully.",
+        )
+
+    def test_unfix_variable(self, simple_supply):
+        simple_supply.backend.fix_variable("energy_cap")
+        simple_supply.backend.unfix_variable("energy_cap")
+        fixed = simple_supply.backend._apply_func(
+            self._is_fixed, simple_supply.backend.variables.energy_cap
+        )
+        assert not fixed.where(fixed.notnull()).all()
+
+    def test_unfix_variable_where(self, simple_supply):
+        where = simple_supply.inputs.energy_cap_max.notnull()
+        simple_supply.backend.fix_variable("energy_cap")
+        simple_supply.backend.unfix_variable("energy_cap", where=where)
+        fixed = simple_supply.backend._apply_func(
+            self._is_fixed, simple_supply.backend.variables.energy_cap
+        )
+        simple_supply.backend.unfix_variable("energy_cap")  # reset
+        assert fixed.sel(techs="test_demand_elec").all()
+        assert not fixed.where(where).all()
diff --git a/tests/test_core_model.py b/tests/test_core_model.py
index 1f011d692..01d8a042e 100644
--- a/tests/test_core_model.py
+++ b/tests/test_core_model.py
@@ -261,46 +261,6 @@ def test_run_mode_mismatch(self):
         )
 
 
-class TestVerboseStrings:
-    def test_verbose_strings_not_implemented(self):
-        m = build_model(
-            {},
-            "simple_supply,two_hours,investment_costs",
-        )
-        with pytest.raises(NotImplementedError) as excinfo:
-            m.verbose_strings()
-
-        assert check_error_or_warning(
-            excinfo,
-            "Call `build()` to generate an optimisation problem before calling this function.",
-        )
-
-    def test_verbose_strings(self, simple_supply):
-        def _compare_to_string(group, component, dims, verbose):
-            component_obj = simple_supply.backend._dataset[component]
-            if verbose:
-                dim_list = ", ".join(dims[k] for k in component_obj.dims)
-                expected = f"{group}[{component}][{dim_list}]"
-            else:
-                expected = f"{group}[{component}][0]"
-            return component_obj.sel(**dims).item().to_string() == expected
-
-        dims_param = {
-            "nodes": "a",
-            "techs": "test_demand_elec",
-            "timesteps": "2005-01-01 00:00",
-        }
-        dims_var = {"carriers": "electricity", **dims_param}
-
-        assert _compare_to_string("parameters", "resource", dims_param, False)
-        assert _compare_to_string("variables", "carrier_con", dims_var, False)
-
-        simple_supply.verbose_strings()
-
-        assert _compare_to_string("parameters", "resource", dims_param, True)
-        assert _compare_to_string("variables", "carrier_con", dims_var, True)
-
-
 class TestValidateMathDict:
     def test_base_math(self, caplog, simple_supply):
         with caplog.at_level(logging.INFO, logger=LOGGER):
diff --git a/tests/test_core_util.py b/tests/test_core_util.py
index 57fe8420e..82651dc31 100644
--- a/tests/test_core_util.py
+++ b/tests/test_core_util.py
@@ -9,13 +9,7 @@
 import pytest
 from calliope.core.util.generate_runs import generate_runs
 from calliope.core.util.logging import log_time
-from calliope.core.util.tools import (
-    copy_docstring,
-    memoize,
-    memoize_instancemethod,
-    validate_dict,
-)
-from calliope.exceptions import ModelError
+from calliope.core.util.tools import memoize, memoize_instancemethod, validate_dict
 
 from .common.util import check_error_or_warning, python36_or_higher
 
@@ -56,19 +50,28 @@ def test_set_log_verbosity(self):
 
         assert logging.getLogger("calliope").getEffectiveLevel() == 50
         assert logging.getLogger("py.warnings").getEffectiveLevel() == 50
-        assert logging.getLogger("calliope.backend.backends").getEffectiveLevel() == 10
+        assert (
+            logging.getLogger("calliope.backend.backend_model").getEffectiveLevel()
+            == 10
+        )
 
         calliope.set_log_verbosity("CRITICAL", include_solver_output=False)
 
         assert logging.getLogger("calliope").getEffectiveLevel() == 50
         assert logging.getLogger("py.warnings").getEffectiveLevel() == 50
-        assert logging.getLogger("calliope.backend.backends").getEffectiveLevel() == 50
+        assert (
+            logging.getLogger("calliope.backend.backend_model").getEffectiveLevel()
+            == 50
+        )
 
         calliope.set_log_verbosity()
 
         assert logging.getLogger("calliope").getEffectiveLevel() == 20
         assert logging.getLogger("py.warnings").getEffectiveLevel() == 20
-        assert logging.getLogger("calliope.backend.backends").getEffectiveLevel() == 10
+        assert (
+            logging.getLogger("calliope.backend.backend_model").getEffectiveLevel()
+            == 10
+        )
 
     def test_timing_log(self):
         timings = {"model_creation": datetime.datetime.now()}
@@ -139,20 +142,6 @@ def test_data_variables_can_be_exported_to_pandas(self, model, variable_name):
             pass
 
 
-class TestCopyDocstring:
-    def test_copy_docstring(self):
-        def _func_w_docstring():
-            "foobar"
-            pass
-
-        def _func(foo, bar):
-            return foo + bar
-
-        docified_func = copy_docstring(_func_w_docstring)(_func)
-        assert docified_func.__doc__ == "foobar"
-        assert docified_func(1, 2) == 3
-
-
 class TestValidateDict:
     @pytest.mark.parametrize(
         ["schema", "expected_path"],
@@ -195,7 +184,7 @@ def test_invalid_dict(self, to_validate, expected_path):
             },
             "additionalProperties": False,
         }
-        with pytest.raises(ModelError) as err:
+        with pytest.raises(calliope.exceptions.ModelError) as err:
             validate_dict(to_validate, schema, "foobar")
         assert check_error_or_warning(
             err,
diff --git a/tests/test_io.py b/tests/test_io.py
index 13b453c6c..ed6ae637d 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -182,7 +182,7 @@ def test_save_read_solve_save_netcdf(self, model, tmpdir_factory):
     def test_save_lp(self, model):
         with tempfile.TemporaryDirectory() as tempdir:
             out_path = os.path.join(tempdir, "model.lp")
-            model.to_lp(out_path)
+            model.backend.to_lp(out_path)
 
             with open(out_path, "r") as f:
                 assert "variables(energy_cap)" in f.read()