Nested constraints and expressions

[brynpickering]

What can be improved?

There are currently two approaches to defining constraints / expressions in the spec:

1. One expression + mask combo for a named expression/constraint

   constraints
     storage_balance_first_ts:
       foreach: [name, snapshot]
       mask: not cyclic_state_of_charge==True and snapshot == get_val_at_index(snapshot=0)
       expression: $soc + ($dispatch * snapshot_weightings) = state_of_charge_initial
     all_other_ts:     
        foreach: [name, snapshot]
        mask: (not cyclic_state_of_charge and not snapshot == get_val_at_index(snapshot=0)) or cyclic_state_of_charge
        expression: $soc + ($dispatch * snapshot_weightings) = roll($soc * (standing_loss ** snapshot_weightings), snapshot=1)

2. Possibility of a list of expression + mask combos, where the masks must be disjoint

   constraints
     storage_balance:
       foreach: [name, snapshot]
       equations: 
       - mask: not cyclic_state_of_charge==True and snapshot == get_val_at_index(snapshot=0)
         expression: $soc + ($dispatch * snapshot_weightings) = state_of_charge_initial
        - mask: (not cyclic_state_of_charge and not snapshot == get_val_at_index(snapshot=0)) or cyclic_state_of_charge
          expression: $soc + ($dispatch * snapshot_weightings) = roll($soc * (standing_loss ** snapshot_weightings), snapshot=1)

(2) is where we got to with the original spec but @FBumann has tended more towards (1).

The reason for (2) was to keep the number of resulting named constraints/expressions small. There are lots of cases where you have a general concept (e.g., balance storage device levels in each timestep) that comes with several variants (what to do in the first timestep differently to all other timesteps). Rather than ask the user to come up with a name for each, they can just provide lists of each variant. It also leads to much more dense expression/constraint arrays which, when storing the definition in an xarray dataset, is beneficial to reduce memory footprint of the built optimisation problem.

However, the list-of-variants approach has a few downsides:

1. it's confusing, hard to debug, and not super readable nor easy to document
2. overriding by merging YAMLs requires overriding all lists at once as you don't have a named key to check against

One downside both approaches have is that you will always want to define variants with disjoint masks and we have no way of knowing that they're disjoint until we evaluate. It's not something we can catch when parsing the spec.

A middle ground might be to enforce naming of variants. As with current list-of-variants approach, these would be rendered as cases of constraints in the latex math and would live in the same, single expression array of expressions.

constraints
  storage_balance:
    foreach: [name, snapshot]
    cases:
      first_ts:
        mask: not cyclic_state_of_charge==True and snapshot == get_val_at_index(snapshot=0)
        expression: $soc + ($dispatch * snapshot_weightings) = state_of_charge_initial
      all_other_ts:     
        mask: (not cyclic_state_of_charge and not snapshot == get_val_at_index(snapshot=0)) or cyclic_state_of_charge
        expression: >-
          $soc + ($dispatch * snapshot_weightings) =
          roll($soc * (standing_loss ** snapshot_weightings), snapshot=1)

constraint groups

constraint_groups:
  storage_balance: [ ... ]

Version

v0.0.0

[FBumann]

This is sth I thought about too, and have a grounded opinion why (1) is the better, simpler approach.

My reasons:

1. The yaml should be ad readable and discussing as possible. each expression having its own name makes this easier, not harder.
2. Accessing duals must be easy. The ID of nested constraints is not well defined
3. The yaml (or dict) should be easy to handle and manipulate. This is much harder with lists.

Accepted Downsides:

1. The user needs to name each constraint

Your proposed middle ground undfortunately does keep the downside, but still makes for an ambiguous dual ID.

[brynpickering]

I think the "variants" approach actually works better for accessing duals since a constraint generally follows comparable logic. If you have three variants for storage dispatch limits then it would be easier to request storage_dispatch_balance.dual than storage_balance_first_ts.dual | storage_balance_middle_ts.dual | storage_balance_final_ts.dual.

[FBumann]

The access seems nicer. But its not a single constraint anymore, is it? So what does this return?

[FBumann]

@brynpickering After some thought I'm open to discuss how to implement constraint groups.

However, should discuss more about in which cases grouping is actually a capability, and improves the model.

I would argue that if we can statically prove that a constraint divided by where/logic atoms covers each row of a constraint exactly once, thats where groups are both mathematically well defined and can be converted to latex, and also are ONE CONSTRAINT.

Your example is statically provable to cover each row exactly once without bowing any data.

How the syntax to express this should be the afterthought.

[FabianHofmann]

before deciding on this in #33 I would like to make another proposal, following the syntax:

expressions:
  previous_soc:
    description: the level a unit carries into a snapshot, by the regime it is in
    branches:
      - where: "cyclic"
        expression: shift(soc, over=snapshot, offset=1, edge='wrap')
      - where: "NOT cyclic AND snapshot > index(snapshot, 0)"
        expression: shift(soc, over=snapshot, offset=1)
      - where: "NOT cyclic AND snapshot == index(snapshot, 0)"
        expression: soc_initial

constraints:
  Storage_energy_balance:
    foreach: [snapshot, storage]
    expression: soc == previous_soc + p_store - p_dispatch

The branches as I call them here (not to confuse things) run on expresssions instead of constraints. this make the constraint itself very readable and underpins the motivation of this original requirement which is branches of one constraint represent the same nature of constraint, here soc == previous_soc + p_store - p_dispatch.
Secondly, you can think of examples where the branches of two added expressions are different: expression a has 2 branches across the time axis, the expression b has three branches across the name axis. The upper proposal would still read simple with the expression masks above leading to 5 where statements. Having the cases in the constraint itself would need 6 where statements. Here is a more complex example from pypsa's define_ramp_limit_constraints . There are three independent axes : the horizon boundary, is_modular, and is_committable.

expressions:
  previous_p:
    branches:
      boundary: {where: "snapshot == index(snapshot, 0)", expression: p_init}
      interior: {where: "snapshot > index(snapshot, 0)", expression: shift(p, over=snapshot, offset=1)}

  ramp_capacity:
    description: what the limit is measured against — a module size where the unit is modular
    branches:
      modular: {where: "is_modular", expression: p_nom_mod}
      whole:   {where: "NOT is_modular", expression: p_nom}

  previous_status:
    branches:
      always_on: {where: "NOT is_committable", expression: 1}
      boundary:  {where: "is_committable AND snapshot == index(snapshot, 0)", expression: status_initial}
      interior:  {where: "is_committable AND snapshot > index(snapshot, 0)", expression: shift(status, over=snapshot, offset=1)}

constraints:
  Generator_ramp_limit_up:
    foreach: [snapshot, generator]
    expression: p - previous_p <= ramp_limit_up * ramp_capacity * previous_status

vs.

constraints:
  Generator_ramp_limit_up:
    foreach: [snapshot, generator]
    cases:
      - where: "snapshot > index(snapshot, 0) AND NOT is_modular AND NOT is_committable"
        expression: p - shift(p, over=snapshot, offset=1) <= ramp_limit_up * p_nom

      - where: "snapshot > index(snapshot, 0) AND NOT is_modular AND is_committable"
        expression: >-
          p - shift(p, over=snapshot, offset=1)
          <= ramp_limit_up * p_nom * shift(status, over=snapshot, offset=1)
          + ramp_limit_start * p_nom * (status - shift(status, over=snapshot, offset=1))

      - where: "snapshot > index(snapshot, 0) AND is_modular AND NOT is_committable"
        expression: p - shift(p, over=snapshot, offset=1) <= ramp_limit_up * p_nom_mod

      - where: "snapshot > index(snapshot, 0) AND is_modular AND is_committable"
        expression: >-
          p - shift(p, over=snapshot, offset=1)
          <= ramp_limit_up * p_nom_mod * shift(status, over=snapshot, offset=1)
          + ramp_limit_start * p_nom_mod * (status - shift(status, over=snapshot, offset=1))

      - where: "snapshot == index(snapshot, 0) AND NOT is_modular AND NOT is_committable"
        expression: p - p_init <= ramp_limit_up * p_nom

      - where: "snapshot == index(snapshot, 0) AND NOT is_modular AND is_committable"
        expression: >-
          p - p_init <= ramp_limit_up * p_nom * status_initial
          + ramp_limit_start * p_nom * (status - status_initial)

      - where: "snapshot == index(snapshot, 0) AND is_modular AND NOT is_committable"
        expression: p - p_init <= ramp_limit_up * p_nom_mod

      - where: "snapshot == index(snapshot, 0) AND is_modular AND is_committable"
        expression: >-
          p - p_init <= ramp_limit_up * p_nom_mod * status_initial
          + ramp_limit_start * p_nom_mod * (status - status_initial)

👀

[FBumann]

@FabianHofmann I like the concept of having the concept on expressions, as it's strictly more flexible.

I would say we should keep where reserved for the top level filter, and rather use another word like when
It's good to make it clear that it's not creating absence like where does.
And it pairs nicely with "cases".

Furthermore, I would probably scope this purely to expressions for now.
And add a mandatory foreach to expressions so proper partition can be proven.
And leave where as a row filter only on constraints, not on expressions

expressions:
  previous_status:
    description: the commitment state a unit carries into a snapshot
    foreach: [snapshot, generator]
    cases:
      always_on: { when: "NOT is_committable", expression: 1 }
      boundary: { when: "is_committable AND position(snapshot) == 0", expression: status_initial }
      interior: { when: "is_committable AND position(snapshot) > 0", expression: shift(status, over=snapshot, offset=1) }

[FabianHofmann]

@FBumann great, good point about the foreach!