Defining piecewise linear constraints¶

In this tutorial, we use the national scale example model to implement a piecewise linear constraint. This constraint will represent a non-linear relationship between capacity and cost per unit capacity of Concentrating Solar Power (CSP).

In [1]:

import numpy as np
import plotly.express as px

import calliope

calliope.set_log_verbosity("INFO", include_solver_output=False)

import numpy as np import plotly.express as px import calliope calliope.set_log_verbosity("INFO", include_solver_output=False)

Model setup¶

Defining our piecewise curve¶

In the base national scale model, the CSP has a maximum rated capacity of 10,000 kW and a cost to invest in that capacity of 1000 USD / kW.

In our updated model, the cost to invest in capacity will vary from 5000 USD / kW to 500 USD / kW as the CSP capacity increases:

In [2]:

capacity_steps = [0, 2500, 5000, 7500, 10000]
cost_steps = [0, 3.75e6, 6e6, 7.5e6, 8e6]

cost_per_cap = np.nan_to_num(np.divide(cost_steps, capacity_steps)).astype(int)

fig = px.line(
    x=capacity_steps,
    y=cost_steps,
    labels={"x": "Capacity (kW)", "y": "Investment cost (USD)"},
    markers="o",
    range_y=[0, 10e6],
    text=[f"{i} USD/kW" for i in cost_per_cap],
)
fig.update_traces(textposition="top center")
fig.show()

capacity_steps = [0, 2500, 5000, 7500, 10000] cost_steps = [0, 3.75e6, 6e6, 7.5e6, 8e6] cost_per_cap = np.nan_to_num(np.divide(cost_steps, capacity_steps)).astype(int) fig = px.line( x=capacity_steps, y=cost_steps, labels={"x": "Capacity (kW)", "y": "Investment cost (USD)"}, markers="o", range_y=[0, 10e6], text=[f"{i} USD/kW" for i in cost_per_cap], ) fig.update_traces(textposition="top center") fig.show()

[2024-12-16 11:51:26] WARNING  /tmp/ipykernel_2989/173276940.py:4: RuntimeWarning: invalid value encountered in divide
  cost_per_cap = np.nan_to_num(np.divide(cost_steps, capacity_steps)).astype(int)

We can then provide this data when we load our model:

Note

We must index our piecewise data over "breakpoints".

In [3]:

new_params = f"""
  parameters:
    capacity_steps:
      data: {capacity_steps}
      index: [0, 1, 2, 3, 4]
      dims: "breakpoints"
    cost_steps:
      data: {cost_steps}
      index: [0, 1, 2, 3, 4]
      dims: "breakpoints"
"""
print(new_params)
new_params_as_dict = calliope.AttrDict.from_yaml_string(new_params)
m = calliope.examples.national_scale(override_dict=new_params_as_dict)

new_params = f""" parameters: capacity_steps: data: {capacity_steps} index: [0, 1, 2, 3, 4] dims: "breakpoints" cost_steps: data: {cost_steps} index: [0, 1, 2, 3, 4] dims: "breakpoints" """ print(new_params) new_params_as_dict = calliope.AttrDict.from_yaml_string(new_params) m = calliope.examples.national_scale(override_dict=new_params_as_dict)

  parameters:
    capacity_steps:
      data: [0, 2500, 5000, 7500, 10000]
      index: [0, 1, 2, 3, 4]
      dims: "breakpoints"
    cost_steps:
      data: [0, 3750000.0, 6000000.0, 7500000.0, 8000000.0]
      index: [0, 1, 2, 3, 4]
      dims: "breakpoints"

[2024-12-16 11:51:26] INFO     Model: initialising

[2024-12-16 11:51:26] INFO     Model: preprocessing stage 1 (model_run)

[2024-12-16 11:51:29] INFO     Model: preprocessing stage 2 (model_data)

[2024-12-16 11:51:29] INFO     Model: preprocessing complete

In [4]:

m.inputs.capacity_steps

m.inputs.capacity_steps

Out[4]:

<xarray.DataArray 'capacity_steps' (breakpoints: 5)> Size: 40B
array([    0.,  2500.,  5000.,  7500., 10000.])
Coordinates:
  * breakpoints  (breakpoints) int64 40B 0 1 2 3 4
Attributes:
    is_result:  False

xarray.DataArray

'capacity_steps'

• breakpoints: 5

• 0.0 2.5e+03 5e+03 7.5e+03 1e+04

  array([    0.,  2500.,  5000.,  7500., 10000.])

• Coordinates: (1)
  • breakpoints
    (breakpoints)
    int64
    0 1 2 3 4

    array([0, 1, 2, 3, 4])

• Indexes: (1)
  • breakpoints
    PandasIndex

    PandasIndex(Index([0, 1, 2, 3, 4], dtype='int64', name='breakpoints'))

• Attributes: (1)

  is_result :

  False

In [5]:

m.inputs.cost_steps

m.inputs.cost_steps

Out[5]:

<xarray.DataArray 'cost_steps' (breakpoints: 5)> Size: 40B
array([      0., 3750000., 6000000., 7500000., 8000000.])
Coordinates:
  * breakpoints  (breakpoints) int64 40B 0 1 2 3 4
Attributes:
    is_result:  False

xarray.DataArray

'cost_steps'

• breakpoints: 5

• 0.0 3.75e+06 6e+06 7.5e+06 8e+06

  array([      0., 3750000., 6000000., 7500000., 8000000.])

• Coordinates: (1)
  • breakpoints
    (breakpoints)
    int64
    0 1 2 3 4

    array([0, 1, 2, 3, 4])

• Indexes: (1)
  • breakpoints
    PandasIndex

    PandasIndex(Index([0, 1, 2, 3, 4], dtype='int64', name='breakpoints'))

• Attributes: (1)

  is_result :

  False

Creating our piecewise constraint¶

We create the piecewise constraint by linking decision variables to the piecewise curve we have created. In this example, we need:

1. a new decision variable for investment costs that can take on the value defined by the curve at a given value of flow_cap;
2. to link that decision variable to our total cost calculation; and
3. to define the piecewise constraint.

In [6]:

new_math = """
  variables:
    piecewise_cost_investment:
      description: "Investment cost that increases monotonically"
      foreach: ["nodes", "techs", "carriers", "costs"]
      where: "[csp] in techs"
      bounds:
        min: 0
        max: .inf
      default: 0
  global_expressions:
    cost_investment_flow_cap:
      equations:
        - expression: "$cost_sum * flow_cap"
          where: "NOT [csp] in techs"
        - expression: "piecewise_cost_investment"
          where: "[csp] in techs"
  piecewise_constraints:
    csp_piecewise_costs:
      description: "Set investment costs values along a piecewise curve using special ordered sets of type 2 (SOS2)."
      foreach: ["nodes", "techs", "carriers", "costs"]
      where: "piecewise_cost_investment"
      x_expression: "flow_cap"
      x_values: "capacity_steps"
      y_expression: "piecewise_cost_investment"
      y_values: "cost_steps"
"""

new_math = """ variables: piecewise_cost_investment: description: "Investment cost that increases monotonically" foreach: ["nodes", "techs", "carriers", "costs"] where: "[csp] in techs" bounds: min: 0 max: .inf default: 0 global_expressions: cost_investment_flow_cap: equations: - expression: "$cost_sum * flow_cap" where: "NOT [csp] in techs" - expression: "piecewise_cost_investment" where: "[csp] in techs" piecewise_constraints: csp_piecewise_costs: description: "Set investment costs values along a piecewise curve using special ordered sets of type 2 (SOS2)." foreach: ["nodes", "techs", "carriers", "costs"] where: "piecewise_cost_investment" x_expression: "flow_cap" x_values: "capacity_steps" y_expression: "piecewise_cost_investment" y_values: "cost_steps" """

Building and checking the optimisation problem¶

With our piecewise constraint defined, we can build our optimisation problem and inject this new math.

In [7]:

new_math_as_dict = calliope.AttrDict.from_yaml_string(new_math)
m.build(add_math_dict=new_math_as_dict)

new_math_as_dict = calliope.AttrDict.from_yaml_string(new_math) m.build(add_math_dict=new_math_as_dict)

[2024-12-16 11:51:29] INFO     Model: backend build starting

[2024-12-16 11:51:29] INFO     Math preprocessing | added file 'plan'.

[2024-12-16 11:51:29] INFO     Math preprocessing | validated math against schema.

[2024-12-16 11:51:29] INFO     Optimisation Model | parameters | Generated.

[2024-12-16 11:51:32] INFO     Optimisation Model | Validated math strings.

[2024-12-16 11:51:32] INFO     Optimisation Model | variables | Generated.

[2024-12-16 11:51:33] INFO     Optimisation Model | global_expressions | Generated.

[2024-12-16 11:51:35] INFO     Optimisation Model | constraints | Generated.

[2024-12-16 11:51:35] INFO     Optimisation Model | piecewise_constraints | Generated.

[2024-12-16 11:51:36] INFO     Optimisation Model | objectives | Generated.

[2024-12-16 11:51:36] INFO     Model: backend build complete

And we can see that our piecewise constraint exists in the built optimisation problem "backend"

In [8]:

m.backend.verbose_strings()
m.backend.get_piecewise_constraint("csp_piecewise_costs").to_series().dropna()

m.backend.verbose_strings() m.backend.get_piecewise_constraint("csp_piecewise_costs").to_series().dropna()

Out[8]:

nodes      techs  carriers  costs   
region1_1  csp    power     monetary    piecewise_constraints[csp_piecewise_costs][reg...
region1_2  csp    power     monetary    piecewise_constraints[csp_piecewise_costs][reg...
region1_3  csp    power     monetary    piecewise_constraints[csp_piecewise_costs][reg...
Name: csp_piecewise_costs, dtype: object

Solve the optimisation problem¶

Once we have all of our optimisation problem components set up as we desire, we can solve the problem.

In [9]:

m.solve()

m.solve()

[2024-12-16 11:51:36] INFO     Optimisation model | starting model in plan mode.

[2024-12-16 11:51:37] INFO     Backend: solver finished running. Time since start of solving optimisation problem: 0:00:01.248357

[2024-12-16 11:51:37] WARNING  Postprocessing: All values < 1e-10 set to 0 in flow_out, flow_in, storage, flow_out_inc_eff, flow_in_inc_eff, cost_operation_variable, cost

[2024-12-16 11:51:37] INFO     Postprocessing: ended. Time since start of solving optimisation problem: 0:00:01.380807

[2024-12-16 11:51:37] INFO     Backend: model solve completed. Time since start of solving optimisation problem: 0:00:01.383840

The results are stored in m._model_data and can be accessed by the public property m.results

Analysing the outputs¶

In [10]:

# Absolute
csp_cost = m.results.cost_investment_flow_cap.sel(techs="csp")
csp_cost.to_series().dropna()

# Absolute csp_cost = m.results.cost_investment_flow_cap.sel(techs="csp") csp_cost.to_series().dropna()

Out[10]:

nodes      carriers  costs   
region1_1  power     monetary    8000000.0
region1_2  power     monetary          0.0
region1_3  power     monetary    3780385.9
Name: cost_investment_flow_cap, dtype: float64

In [11]:

# Relative to capacity
csp_cap = m.results.flow_cap.sel(techs="csp")
csp_cost_rel = csp_cost / csp_cap
csp_cost_rel.to_series().dropna()

# Relative to capacity csp_cap = m.results.flow_cap.sel(techs="csp") csp_cost_rel = csp_cost / csp_cap csp_cost_rel.to_series().dropna()

Out[11]:

nodes      carriers  costs   
region1_1  power     monetary     800.00000
region1_3  power     monetary    1492.00507
dtype: float64

In [12]:

# Plotted on our piecewise curve
fig.add_scatter(
    x=csp_cap.to_series().dropna().values,
    y=csp_cost.to_series().dropna().values,
    mode="markers",
    marker_symbol="cross",
    marker_size=10,
    marker_color="cyan",
    name="Installed capacity",
)
fig.show()

# Plotted on our piecewise curve fig.add_scatter( x=csp_cap.to_series().dropna().values, y=csp_cost.to_series().dropna().values, mode="markers", marker_symbol="cross", marker_size=10, marker_color="cyan", name="Installed capacity", ) fig.show()

Troubleshooting¶

If you are failing to load a piecewise constraint or it isn't working as expected, here are some common things to note:

1. The extent of your x_values and y_values will dictate the maximum values of your piecewise decision variables. In this example, we define capacity_steps over the full capacity range that we allow our CSP to cover in the model. However, if we set capacity_steps to [0, 2500, 5000, 7500, 9000] then flow_cap would never go above a value of 9000.

2. The x_values and y_values parameters must have the same number of breakpoints and be indexed over breakpoints. It is possible to extend these parameters to be indexed over other dimensions (e.g., different technologies with different piecewise curves) but it must always include the breakpoints dimension.

3. x_values must increase monotonically. That is, [0, 5000, 2500, 7500, 10000] is not valid for capacity_steps in this example. y_values, on the other hand, can vary any way you like; [0, 6e6, 3.75e6, 8e6, 7.5e6] is valid for cost_steps.

4. x_expression and y_expression must include reference to at least one decision variable. It can be a math expression, not only a single decision variable. flow_cap + storage_cap / 2 would be valid for x_expression in this example.

5. Piecewise constraints will make your problem more difficult to solve since each breakpoint adds a binary decision variable. Larger models with detailed piecewise constraints may not solve in a reasonable amount of time.