Special Ordered Sets of type 2 (SOS2) for piecewise linearisation

[brynpickering]

With the new math formulation, it is possible to define simple piecewise linear constraints (e.g. for investment costs), possibly solving #107.

In these, a convex hull is created that approximates the non-linear curve and the structure of the model ensures that the variable will sit somewhere along the upper bound of that hull and therefore take on a value very similar to the nonlinear curve's values (see image, Fig 4.2 from my phd).

For example, efficiency might change with technology output. As the output increases on the Y axis, the curve lead to increasingly more input requirements, therefore simulating a reduced efficiency with greater input. This works because the model wants to maximise output for each unit input (otherwise it has to "spend" more on the input carrier, increasing the objective function value) and so will have to sit on the upper bounding curve.

The advantage of this approach is that the model can avoid lots of integer/binary variables. The disadvantage is that it can only be used to describe problems where things get worse with increased output/capacity of a technology. In reality, economies of scale lead to lower investment costs for increased capacity, and part-load efficiencies are often lower than full load efficiency.

To describe these problems, I believe we need to use Special Ordered Sets. in SOS of type 2 - SOS2 - each "joint" on the piecewise curve is described by a binary variable and only two adjacent variables may be non-zero at any time. This means only one section of a piecewise curve is "active" at a given moment and so input and output must be related according to the slope of the line between those two points. SOS2 come with the advantage of being able to describe most non-linear curve shapes in the form of a discrete set of straight lines. But they need many more binary variables so might make most problems intractable. Nevertheless, this should be an option for Calliope modellers.

Solvers have python APIs to define an SOS that we could tap into (Pyomo, Gurobi). They also have more generalised (and easier to use) piecewise linear constraint interfaces (Pyomo, Gurobi) but we still need to provide a way for users to construct them.

If we use the generalised piecewise linear constraint constructor, I would force it to use SOS2 with the "==" operator between x and y. This is partly because it is all that Gurobi can handle (Pyomo has many more options), but also because it is the most likely formulation of a piecewise constraint that someone will think of.

We can't use the standard math syntax we have. Instead we need something that allows reference to the x and y variables being compared and the values of x and y at each "breakpoint". Some ideas:

1

piecewise_constraints:
  my_constraint:
    foreach: [nodes, techs]
    x: 
       variable: flow_cap
       values: [0, 1, 2, 3, 4]
    y: 
       variable: piecewise_cost_investment
       values: [0, 10, 15, 18, 20]

2

piecewise_constraints:
  my_constraint:
    foreach: [nodes, techs]
    x: flow_cap
    x_values: [0, 1, 2, 3, 4]
    y: piecewise_cost_investment
    y_values: [0, 10, 15, 18, 20]

3

piecewise_constraints:
  my_constraint:
    foreach: [nodes, techs]
    x: flow_cap
    y: piecewise_cost_investment
    breakpoints: 
      0: 0
      1: 10
      2: 15
      3: 18
      4: 20

4

piecewise_constraints:
  my_constraint:
    foreach: [nodes, techs]
    x: flow_cap
    y: piecewise_cost_investment
    breakpoints: 
      - [0, 0]
      - [1, 10]
      - [2, 15]
      - [3, 18]
      - [4, 20]

5

piecewise_constraints:
  my_constraint:
    foreach: [nodes, techs]
    x: flow_cap
    y: piecewise_cost_investment
    breakpoints:       
      - {x: 0, y: 0}
      - {x: 1, y: 10}
      - {x: 2, y: 15}
      - {x: 3, y: 18}
      - {x: 4, y: 20}

Option 1 and 2 are probably better as breakpoints will likely be defined per technology / node and that could be defined by an input parameter in those options. E.g.,:

1:

piecewise_constraints:
  my_constraint:
    foreach: [nodes, techs]
    x: 
       variable: flow_cap
       values: my_constraint_breakpoint_data_x
    y: 
       variable: piecewise_cost_investment
       values: my_constraint_breakpoint_data_y

and...

techs:
  some_tech:
    my_constraint_breakpoint_data_x: 
      data: [0, 1, 2, 3, 4]
      index: [0, 1, 2, 3, 4]
      dims: breakpoints
    my_constraint_breakpoint_data_y:
      data: [0, 1, 2, 3, 4]
      index: [0, 10, 15, 18, 20]
      dims: breakpoints

However, there would need to be an additional model dimension for breakpoints that this parameter is indexed over. All the options would work if you have to define one piecewise constraint for each index item in foreach. Maybe that's OK because it's unlikely that lots of different piecewise constraints are going to be needed (one per tech for 5-15 techs?) but it is a break from the current math format.

Or, there could be one parameter per breakpoint, so the options would become, e.g.:

1:

piecewise_constraints:
  my_constraint:
    foreach: [nodes, techs]
    x: 
       variable: flow_cap
       values: [my_constraint_breakpoint_data_1_x, my_constraint_breakpoint_data_2_x, ...]
    y: 
       variable: piecewise_cost_investment
       values: [my_constraint_breakpoint_data_1_y, my_constraint_breakpoint_data_2_y, ...]

4:

piecewise_constraints:
  my_constraint:
    foreach: [nodes, techs]
    x: flow_cap
    y: piecewise_cost_investment
    breakpoints:
      - [my_constraint_breakpoint_data_1_x, my_constraint_breakpoint_data_1_y]
      - [my_constraint_breakpoint_data_2_x, my_constraint_breakpoint_data_2_y]

[jmorrisnrel]

Has there been any consideration of perhaps farther down the line looking for areas where certain piecewise linear parameters could be linearized by splitting out the nodes/techs into separate linear technologies? A specific example would be increasing LCOE-based cost-capacity curves where economies of scale are outweighed by increasing LCOEs as less-desirable locations are included (assuming there is no base cost/y offset incurred). This would be highly dependent on other constraints fitting around splitting nodes within the model (i.e. a non-linear efficiency curve wouldn't work with a split technology without other processing done/a more complex constraint added), but many of our technical resource assessments output linearized costs for a number of points that would need to be aggregated into one region/node for a Calliope model, with build decisions reliably going from lowest to highest LCOE.

A modeler could preprocess that data to create linear sub-technologies to represent the different pieces of the curve, but if there are opportunities to do some of that within Calliope and increase model size/performance and ease of modeling that would be valuable to some users.

[brynpickering]

LCOE is an output of the model as it is the (total investment in capacity + the total operating costs ) / (total electricity output) that defines it. My interpretation of your case is that the closer you get to the use of all e.g., available area, the less desirable the remaining area is. This might be because the cost of land increases, access for labour and materials is restricted, or because the yield from these locations will be lower. All would adversely affect the LCOE as you have depicted in that curve.

Your figure looks like the combinatorial curve of economies of scale (the flattening off of the curve from zero to 20GW) and then pushing capacity into less desirable locations (the rising curve from 130GW). This could be depicted in Calliope using piecewise linearisation, but would require a different Y-axis. E.g., total investment cost. Perhaps two separate curves would be possible, one for cumulative capacity against investment cost and the other for cumulative capacity against maximum power output (i.e., capacity factor).

I'd be more inclined to say that defining different technologies (expensive_productive_pv, cheap_productive_pv, expensive_not_productive_pv, ...) would be better as it gives you the possibility of more fine-grain control and perhaps makes it easier to interpret the results.

Still, your example is a good one for using SOS2 since the curve changes the sign of its gradient. That's something that the convex hull approach can't handle (see case (d) in my image).

Let me know if I've got the wrong end of the stick as to the problem you're describing!

[eeroinkeri]

Piecewise constraints would be great addition to Calliope! I have been trying to simulate the tutorial: https://calliope.readthedocs.io/en/latest/examples/piecewise_constraints/

However, I haven't managed to get it work

I have tried it two ways:

1. directy via jupyter

m.math["piecewise_constraints"]["csp_piecewise_costs"] = {
    "description": "Set investment costs values along a piecewise curve using special ordered sets of type 2 (SOS2).", 
    "foreach": ["nodes", "techs", "carriers", "costs"], 
    "where": "piecewise_cost_investment", 
    "x_expression": "flow_cap", 
    "x_values": "capacity_steps", 
    "y_expression": "piecewise_cost_investment", 
    "y_values": "cost_steps", 
}

It gives error:

KeyError: 'piecewise_constraints'

2. Additional math through override and yaml file

piecewise_constraints:
  csp_piecewise_costs:
    description: Set investment costs values along a piecewise curve using special ordered sets of type 2 (SOS2).
    foreach: [nodes, techs, carriers, costs]
    where: "[csp] in techs"
    x_expression: flow_cap
    x_values: capacity_steps
    y_expression: piecewise_cost_investment
    y_values: cost_steps

Now all new math are succesfully added to the model.math, but new error appears during model.build():

ModelError: Errors during validation of the math dictionary.:
 * Additional properties are not allowed ('piecewise_constraints' was unexpected)

Should this tutorial work as it is, or do I maybe have some misconfiguration?

[brynpickering]

What version of Calliope are you using? You can check with:

import calliope
calliope.__version__

[eeroinkeri]

Hi, it is 0.7.0.dev3

[brynpickering]

OK, update to v0.7.0.dev4, where the feature was added. You should be able to do mamba update calliope to get the newer version

[eeroinkeri]

Thank you, I need to pay more attention to the version and version history..

[Longquan-Li]

Thanks for the function that I found very useful.

However, I met some issue when I tested it using dev4. It seems that the global_expressions defined in the yaml code provided in dev4, as shown below, does not effectively inject the piecewise_cost_investment into the cost_investment_flow_cap,

cost_investment_flow_cap:
equations:
- expression: "$cost_sum * flow_cap"
where: "NOT [pipeline_new, overhead_line] in techs"
- expression: "piecewise_cost_investment"
where: "[pipeline_new, overhead_line] in techs"

After I saved the results in CSV, I did not see the item related to the techs that use the piece_cost_investment in the file cost_investment_flow_cap.csv. Instead, the piecewise_cost_investment is shown in another separate file. However, in this way, it seems that the piecewise_cost_investment is not taken into account in the optimization objective, according to my tests.

[brynpickering]

@Longquan-Li it's difficult to be sure why this might be without more information. I've checked our docs example for piecewise costs and it does pass the costs onto cost_investment_flow_cap without any issue. Are you able to share a MWE?

[brynpickering]

So, it's because the global expression isn't triggering because it isn't getting past the top-level where string. You need to update the math to:

global_expressions:
  cost_investment_flow_cap:
    where: (flow_cap AND (cost_flow_cap OR cost_flow_cap_per_distance)) OR piecewise_cost_investment
    equations:
      - expression: "$cost_sum * flow_cap"
        where: "NOT [pipeline_new, overhead_line] in techs"
      - expression: "piecewise_cost_investment"
        where: "[pipeline_new, overhead_line] in techs"

[brynpickering]

@Longquan-Li I have updated the docs example to reflect this: #779

[Longquan-Li]

Thank you. I tried the code you shared. The problem was solved, but I came across an error when I tried to save the results as csv.

Error in astype, /Users/longquanli/conda/envs/calliope/lib/python3.12/site-packages/xarray/core/duck_array_ops.py:215
could not convert string to float: '147.0 C0'

I am not sure what's the reason behind this.

In the same while, I tried the global expression in another example, as follows,

global_expressions:
cost_investment:
where: >
(cost_investment_flow_cap OR cost_investment_storage_cap
OR cost_investment_source_cap OR cost_investment_area_use
OR cost_investment_purchase OR piecewise_cost_investment)
equations:
- expression: >
sum(default_if_empty(cost_investment_flow_cap, 0), over=carriers) +
default_if_empty(cost_investment_storage_cap, 0) +
default_if_empty(cost_investment_source_cap, 0) +
default_if_empty(cost_investment_area_use, 0) +
default_if_empty(cost_investment_purchase, 0) +
default_if_empty(piecewise_cost_investment, 0)

The error did not occur and the results seem reasonable.

[Longquan-Li]

@brynpickering , the error disappeared after the code being revised to,

global_expressions:
cost_investment_flow_cap:
where: (flow_cap AND (cost_flow_cap OR cost_flow_cap_per_distance)) OR piecewise_cost_investment
equations:
- expression: "$cost_sum * flow_cap"
where: "NOT [pipeline_new, overhead_line] in techs"
- expression: "1.0 * piecewise_cost_investment"
where: "[pipeline_new, overhead_line] in techs"

[Longquan-Li]

MWE.zip

Sure, here it is.

[jmorrisnrel]

Has there been any success using an SOS2 constraint for technology efficiency without it becoming non-linear? It intuitively makes sense that essentially scaling flow_out with itself would result in a non-linear constraint even with SOS2 piecewise formulation but I'm wondering if someone has come up with a creative set of constraints/expressions/decision variables for defining a non-convex efficiency curve. Capturing reduced efficiencies at lower power output would be valuable for operate mode runs where we have the ability to model more detailed constraints. The closest I've managed to get is something that ties the source_use/flow_in to the capacity with the appropriate value for that point on a heat curve but that doesn't actually scale with flow_out so you get a flat input consumption between breakpoints.

Turns out the SOS2 constraint interpolates between breakpoints and I just misread my test results and confused myself. Leaving this up in case someone else is trying to do the same. The below piecewise constraint will work for supplies, should just need to replace source_use * source_eff with flow_in * flow_in_eff for conversion/transmission and update balance_conversion/balance_transmission instead of balance_supply.

tech:
   flow_out_eff_piecewise_x:
      data: [0, 0.01, 0.5, 0.75, 0.85, 0.95, 1]
      index: [0, 1, 2, 3, 4, 5, 6]
      dims: breakpoints
    flow_out_eff_piecewise_y:
      data: [0, 0.02, 0.66, 0.885, 0.8925, 1.1875, 1.25]
      index: [0, 1, 2, 3, 4, 5, 6]
      dims: breakpoints

piecewise_constraints:
  piecewise_efficiency:
    description: >
      Use SOS2 to model piecewise linear output efficiency curves/heat rates.
    foreach: [nodes, techs, timesteps, carriers]
    where: >
      any(flow_out_eff_piecewise_x, over=breakpoints) AND any(flow_out_eff_piecewise_y, over=breakpoints)
    x_expression: flow_out / flow_cap
    x_values: flow_out_eff_piecewise_x
    y_expression: source_use * source_eff / flow_cap
    y_values: flow_out_eff_piecewise_y

# Override balance_supply to not apply to the piecewise tech as we're handling it in the piecewise constraint
constraints:
  balance_supply_no_storage:
    description: >-
      Fix the outflow of a `supply` technology to its consumption of the available source.
    foreach: [nodes, techs, carriers, timesteps]
    where: "base_tech=supply AND NOT include_storage=True"
    equations:
      - where: NOT flow_out_eff_piecewise
        expression: flow_out_inc_eff == source_use * source_eff

Using any number of options that tie flow_out to the y expression in any way results in a nonlinear term error from Pyomo so you need to explicitly tie source_use/flow_in to the y expression.

piecewise_constraints:
  piecewise_efficiency:
    description: >
      Use SOS2 to model piecewise linear output efficiency curves/heat rates.
    foreach: [nodes, techs, timesteps, carriers]
    where: >
      any(flow_out_eff_piecewise_x, over=breakpoints) AND any(flow_out_eff_piecewise_y, over=breakpoints)
    x_expression: flow_out / flow_cap
    x_values: flow_out_eff_piecewise_x
    y_expression: flow_out / source_use * source_eff
    y_values: flow_out_eff_piecewise_y

ValueError: Model constraint (piecewise_constraints[piecewise_efficiency][0].c[1]) contains nonlinear terms that cannot be written to LP format