Custom constraints

[brynpickering]

In 0.6.0 we have (temporarily?) removed the ability to introduce custom constraints via the YAML format¹.

To bring it back in, ideally a user could have one Python script which introduces any new sets, parameters, decision variables, and constraints. If we provide a template for this, it would ensure that a user includes all the necessary information to produce the desired results.

Steps to introduce the functionality

• [ ] Get list of all data required to build a brand new constraint
• [ ] Build a custom constraint template and test with a new constraint

¹ A user familiar with Pyomo can currently add a custom constraint after building the backend model (model.run(build_only=True)) by accessing the Pyomo ConcreteModel directly via model._backend_model.

[sjpfenninger]

Always sets for each tech and tech_group to allow you to easily attach per-tech custom constraints via those sets?

[timtroendle]

Here's a work-in-progress suggestion of how custom constraints could look like in the yaml file. The idea is that next to the user-friendly Calliope structure, one can manually add Sets, Parameters, Variables, and Constraints.

The constraints defined below are the group constraints in Calliope 0.6.5. It requires #322 to be solved first.

custom_constraints:
    sets:
        wind_techs:
            elements:
                - onshore_wind
                - offshore_wind
            within: techs
        ac_transmission:
            elements:
                - ac_transmission
            within: techs
        nordic_countries:
            elements:
                - FIN
                - CHE
                - NOR
            within: locs
        tims_favorite_timesteps:
            elements:
                - "2016-04-01 04:00"
                - "2016-06-16 08:00"
            within: timesteps
        monetary_cost:
            elements:
                - monetary
            within: costs
    parameters:
        carrier_prod_share: 0.25
        net_import_share: 0.7
        energy_cap_share: 0.2
        energy_cap: 200
        cost_cap: 1e6
        cost_var_cap: 1e4
        cost_investment_cap: 1e3
        resource_area_cap: 100
    # TODO add variables
    constraints:
        net_import_share:
            foreach: [nordic_locs, carriers]
            eq: net_import_sum <= net_import_share * demand_sum
            components:
                net_import_sum:
                    sum: carrier_prod + carrier_con
                    over: [ac_transmission, timesteps]
                demand_sum:
                    sum: carrier_con
                    over: [demand, timesteps]

        carrier_prod_share:
            eq: energy_sum_wind <= carrier_prod_share * total_energy_sum
            foreach: [nordic_locs, carriers]
            components:
                energy_sum_wind:
                    sum: carrier_prod
                    over: [wind_techs, timesteps]
                total_energy_sum:
                    sum: carrier_prod
                    over: [supply_techs, timesteps]
            # TODO: allow for conditional statements

        energy_cap_share:
            eq: energy_cap_sum_wind <= energy_cap_share * energy_cap_sum_all
            foreach: [locs, electricity]
            components:
                energy_cap_sum_wind:
                    sum: energy_cap
                    over: [wind_techs]
                energy_cap_sum_all:
                    sum: energy_cap
                    over: [techs_supply]

        energy_cap:
            eq: energy_cap_sum_wind <= energy_cap
            foreach: [locs, electricity]
            components:
                energy_cap_sum_wind:
                    sum: energy_cap
                    over: [wind_techs]

        carrier_prod_share_per_timestep:
            eq: energy_sum_wind <= carrier_prod_share * total_energy_sum
            foreach: [nordic_locs, carriers, tims_favorite_timesteps]
            components:
                energy_sum_wind:
                    sum: carrier_prod
                    over: [wind_techs]
                total_energy_sum:
                    sum: carrier_prod
                    over: [supply_techs]

        demand_share:
            foreach: [locs, carriers]
            eq: net_import_sum <= demand_share * -1 * demand_sum
            components:
                net_import_sum:
                    sum: carrier_prod
                    over: [wind_techs, timesteps]
                demand_sum:
                    sum: carrier_con
                    over: [demand, timesteps]

        resource_area:
            foreach: []
            eq: resource_area_sum <= resource_area_cap
            components:
                resource_area_sum:
                    sum: resource_area
                    over: [wind_techs, nordic_countries]

        cost_cap:
            foreach: [monetary_cost]
            eq: cost_sum <= cost_cap
            components:
                cost_sum:
                    sum: cost
                    over: [locs, techs, timesteps]

        cost_var_cap:
            foreach: [monetary_cost]
            eq: cost_sum <= cost_cap
            components:
                cost_sum:
                    sum: cost_var
                    over: [locs, techs, timesteps]

        cost_investment_cap:
            foreach: [monetary_cost]
            eq: cost_sum <= cost_cap
            components:
                cost_sum:
                    sum: cost_investment
                    over: [locs, techs]
    #TODO needs a variable
    # c1: carrier_prod[hp, t] / resource[demand_tech, t] == decision_var[hp]
    # c2: decision_var[hp] + decision_var[boiler] = 0.8

[brynpickering]

Some thoughts on this. We can leverage the update to how constraint sets are built in #346 to make it easier to define constraints in YAML. In config/subsets.yaml, each named constraint has the dimensions over which the constraint is built and a where "mask" to tell it whether to consider an indexed component or not. E.g.:

carrier_consumption_max:
        foreach: [nodes, techs, carriers, timesteps]
        where: [carrier, and, [inheritance(transmission), or, inheritance(demand), or, inheritance(storage)], and, [not cap_method='integer', or, inheritance(demand)], and, allowed_carrier_con=True]
        subset:
            carrier_tiers: [in]

This means that the internal constraints should never need to check if the indexed component is valid, since it is pre-filtered when generating the constraint set.

For the above, the YAML constraint should then be able to look like this, without any additional 'components':

carrier_consumption_max:
    eq: carrier_con >= -1 * energy_cap * timestep_resolution

[brynpickering]

For the above, the YAML constraint should then be able to look like this, without any additional 'components':

carrier_consumption_max:
    eq: carrier_con >= -1 * energy_cap * timestep_resolution

On another note, I do wonder if this would be better looking like this:

carrier_consumption_max:
    eq: carrier_con[carrier, node, tech, timestep] >= -1 * energy_cap[node, tech] * timestep_resolution[timestep]

I.e., is the mapping of dimension names to the relevant variables/parameters something we should infer automatically (the first option) or should they be explicitly given?

[brynpickering]

Another thing we've now implemented, that makes this particular addition easier, is the definition of variables. This means that the commented section at the end of @timtroendle's work-in-progress suggestion can be dealt with there.

[brynpickering]

Still, there are some issues I can see coming up:

1. we need a hierarchy of checking constraint components. First, we check if they're an internally defined parameter/variable, second we check if they're a constraint component, third we check that sub-components of constraint components are internally defined parameters/variables.
2. we need to handle conditionals, either by having a way of switching within the constraint definition (sounds messy) or by splitting up constraints first, to avoid any need for conditionals. Splitting up constraints might lead to quite a large number of them (e.g., energy balance constraints have a lot of internal conditionals), but would be something we could handle on the set generation side of things, since that is structured to combine a bunch of conditionals into a valid, constraint-specific subset.
3. If we don't have explicit dimension referencing for the variables/parameters (e.g., carrier_con[carrier, node, tech, timestep]), how do we enable a constraint to be defined that fixes a specific dimension (e.g., we build over all technologies at a specific location)? Perhaps this is a reason to be explicit with dimension referencing, or something to account for in per-constraint components. I.e., as well as sum there could be a way to select a specific dimension value.
4. We might need custom coded logic that doesn't work at all with YAML. One example in line with (3) is that we might want a constraint to focus only on the first timestep of the timeseries (see here for a constraint where that is necessary). However, the first timestep is user-defined (e.g., dependant on the year). So we would need some way of referring to the first element of the set, which would need some special YAML logic OR a simple python script that is referred to.

[brynpickering]

Here are some examples of possible implementations that handle complex core calliope constraints. The idea is that if we can handle these constraints, there is high chance that we can handle most constraints a user might define.

It turns out that (2) above is particularly difficult. Splitting into individual constraints creates a massive amount of repetition of constraint components and doesn't seem tenable. Below are two possible options that could work: using an 'if-else' syntax within the constraint definition and using the 'where' syntax that we introduced in constraint subsetting

1. Balance supply

Using 'if-else' syntax

```yaml balance_supply: foreach: [carrier, node, tech, timestep] where: [resource, and, inheritance(supply)] eq: - if: energy_eff == 0 then: carrier_prod == 0 - if: force_resource == 1 then: carrier_prod_div_energy_eff == available_resource - if: min_use then: min_use * available_resource <= carrier_prod_div_energy_eff <= available_resource - else: then: carrier_prod_div_energy_eff <= available_resource components: available_resource: - if: resource_unit == energy_per_area then: resource * resource_scale * resource_area - if: resource_unit == energy_per_cap then: resource * resource_scale * energy_cap - else: then: resource * resource_scale carrier_prod_div_energy_eff: carrier_prod / energy_eff ```

Using 'where' syntax

```yaml balance_supply: foreach: [node in nodes, tech in techs, carrier in carriers, timestep in timesteps] where: [resource, and, inheritance(supply)] options: - where: [energy_eff=0] eq: carrier_prod[node, tech, carrier, timestep] == 0 - where: [force_resource=True] eq: carrier_prod_div_energy_eff == available_resource - where: [resource_min_use>0] eq: resource_min_use[node, tech] * available_resource <= carrier_prod_div_energy_eff <= available_resource - eq: carrier_prod_div_energy_eff <= available_resource - where: [resource_unit='energy_per_area'] components: available_resource: eq: resource[node, tech, timestep] * resource_scale[node, tech] * resource_area[node, tech] - where: [resource_unit='energy_per_cap'] components: available_resource: eq: resource[node, tech, timestep] * resource_scale[node, tech] * energy_cap[node, tech] - where: [resource_unit='energy'] components: available_resource: eq: resource[node, tech, timestep] * resource_scale[node, tech] components: carrier_prod_div_energy_eff: eq: carrier_prod[node, tech, carrier, timestep] / energy_eff[node, tech, timestep] ```

3. Balance supply plus

Using 'if-else' syntax

```yaml balance_supply_plus: foreach: [node, tech, carrier, timestep] where: [inheritance(supply_plus)] eq: - if: include_storage is False then: available_resource == carrier_prod_incl_losses - if: include_storage is True and storage_inter_cluster and lookup_cluster_first_timestep is True then: storage == available_resource - carrier_prod_incl_losses - else: storage == storage_at_timestep_start + available_resource - carrier_prod_incl_losses components: carrier_prod_incl_losses: eq: - if: energy_eff == 0 or parasitic_eff == 0 then: 0 - else: carrier_prod / (energy_eff * parasitic_eff) available_resource: eq: resource_con * resource_eff storage_at_timestep_start: eq: - if: include_storage=True, and, not run.cyclic_storage=true, and, timestep = timesteps[0] then: storage_initial * storage_cap - else: ((1 - storage_loss) ** timestep_resolution[timesteps=previous_step]) * storage[timesteps=previous_step] index_items: previous_step: eq: - if: include_storage is True and run.cyclic_storage is True and timestep == timesteps[0] then: timesteps[-1] - if: include_storage is True and timestep > timesteps[0] then: timesteps[index(timestep) - 1] # needs a function 'index' to get the position of 'timestep' in the list 'timesteps' <- could be a lookup table, like below - if: include_storage is True and timestep>timesteps[0] and clusters and lookup_cluster_first_timestep is True then: lookup_cluster_last_timestep[timestep] ```

Using 'where' syntax

```yaml balance_supply_plus: foreach: [node in nodes, tech in techs, carrier in carriers, timestep in timesteps] options: - where: [inheritance(supply_plus), and, not include_storage=True] eq: available_resource == carrier_prod_incl_losses - where: [inheritance(supply_plus), and, include_storage=True, and, storage_inter_cluster, and, lookup_cluster_first_timestep=True] eq: storage[node, tech, timestep] == available_resource - carrier_prod_incl_losses - where: [inheritance(supply_plus), and, include_storage=True, and, not run.cyclic_storage=true, and, timestep = timesteps[0]] eq: storage[node, tech, timestep] == storage_initial[node, tech] * storage_cap[node, tech] + available_resource - carrier_prod_incl_losses components: previous_step: elements: timesteps[-1] within: timesteps - where: [inheritance(supply_plus), and, include_storage=True, and, timestep > timesteps[0]] eq: storage[node, tech, timestep] == storage_at_timestep_start + available_resource - carrier_prod_incl_losses components: previous_step: elements: within: timesteps - where: [inheritance(supply_plus), and, include_storage=True, and, timestep > timesteps[0], and, clusters, and, lookup_cluster_first_timestep=True] eq: storage[node, tech, timestep] == storage_at_timestep_start + available_resource - carrier_prod_incl_losses components: previous_step: elements: lookup_cluster_last_timestep[timestep] within: timesteps components: carrier_prod_incl_losses: - where: [energy_eff=0, or, parasitic_eff=0] eq: 0 - where: [energy_eff=0, or, parasitic_eff=0] eq: carrier_prod[node, tech, carrier, timestep] / (energy_eff[node, tech, timestep] * parasitic_eff[node, tech, timestep]) available_resource: eq: resource_con[node, tech, timestep] * resource_eff[node, tech] storage_at_timestep_start: eq: ((1 - storage_loss[node, tech]) ** timestep_resolution[previous_step]) * storage[node, tech, previous_step] ```

[brynpickering]

In the 'balance supply' constraint, storage is there and so is the need for a link between timesteps. One option is to introduce two helper functions (to add to 'inheritance' in subsetting): get_index and get_item. get_index would get the index number of an item in a list. get_item would get the item from a list based on the index position.

[sjpfenninger]

New thinking summarised:

Processing approach

The approach for processing constraints in three steps:

1. Human-readable constraint definitions are parsed into a list of equations, where each list entry is a dict with keys "foreach", "where", "equation"

2. Iterate over the list of equations from (1), safely parsing each equation string from using pyparsing

The final step is backend-specific - here the actual model objects are built:

3. Iterate over the list of parsed equations from (2), building the actual model components with the selected API (i.e., Gurobi, or Pyomo)

Human-readable constraint definition format

The top-level YAML key is "constraints", which is a list of dicts.

Allowed keys:

• foreach
• where (optional)
• one of equation or equations
• components (optional)
• index_items (optional)

The equation can be given as as string (equation) or a list of expression dicts ({"where": ..., "expression": ...}).

The keys equations, components, index_items only accept a list of expression dicts ({"where": ..., "expression": ...}).

For index_items, an additional key is required: on_dimension. This specifies which dimension name in the main expression(s) the index item applies to.

where is always optional: if not given, it is assumed to be an empty list. Thus, if an expression applies without further subsetting, it can either be written as:

- where: []
  expression: ...

or

- expression: ...

[brynpickering]

This is available and documented in v0.7