A horizon too large for one solve¶

A model over a year can be too large to solve at once. Solve it a window at a time instead, and see what that costs. Every block runs from the examples/ folder.

The model¶

tutorial/storage.yaml is the first model with two additions. A generator is either on or off, and gas cannot run below 50 while it is on. A battery carries its level from one snapshot to the next.

tutorial/storage.yaml

description: Least-cost dispatch with on-off generators and a battery.

dimensions:
  snapshot:
    description: dispatch periods
    dtype: int
  generator:
    description: generating units

parameters:
  p_max:
    description: installed capacity
    dims: [generator]
  load:
    description: demand to be met
    dims: [snapshot]
  cost:
    description: marginal cost
    dims: [generator]
  p_min:
    description: least output of a generator that is on
    dims: [generator]
  no_load_cost:
    description: cost of one snapshot on, whatever the output
    dims: [generator]
  soc_initial:
    description: battery level before the first snapshot
    dims: []

variables:
  p:
    description: output of a generator in a snapshot
    dims: [snapshot, generator]
    where: "p_max > 0"
    bounds:
      lower: 0
      upper: p_max
  on:
    description: whether a generator is on in a snapshot
    dims: [snapshot, generator]
    where: "p_max > 0"
    domain: binary
  charge:
    description: energy into the battery
    dims: [snapshot]
    bounds:
      lower: 0
      upper: 30
  discharge:
    description: energy out of the battery
    dims: [snapshot]
    bounds:
      lower: 0
      upper: 30
  soc:
    description: battery level at the end of a snapshot
    dims: [snapshot]
    bounds:
      lower: 0
      upper: 100

constraints:
  power_balance:
    dims: [snapshot]
    expression: sum(p, over=generator) + discharge - charge == load
  p_upper:
    description: a generator that is off produces nothing
    dims: [snapshot, generator]
    expression: p <= p_max * on
  p_lower:
    description: a generator that is on produces at least its minimum
    dims: [snapshot, generator]
    expression: p >= p_min * on
  soc_first:
    description: the first snapshot starts from soc_initial
    dims: [snapshot]
    where: "snapshot == 0"
    expression: soc == soc_initial + 0.9 * charge - discharge
  soc_next:
    description: every other snapshot starts from the level the one before ended with
    dims: [snapshot]
    expression: soc == shift(soc, along=snapshot, offset=1) + 0.9 * charge - discharge

objective:
  sense: minimize
  description: total cost of generation and of keeping generators on
  expression: sum(p * cost) + sum(on * no_load_cost)

Solve one day. The sources gain p_min, no_load_cost and the battery's level at the start, soc_initial:

import polars as pl

import specsolve as sps

generators = pl.read_csv('tutorial/data/generators.csv')
load = pl.read_csv('tutorial/data/load.csv')
fleet = {'generator': generators['generator']} | {
    name: generators.select('generator', value=name) for name in ('p_max', 'cost', 'p_min', 'no_load_cost')
}

result = sps.solve('tutorial/storage.yaml', fleet | {'snapshot': load['snapshot'], 'load': load, 'soc_initial': 0.0})
print(result.objective)

10850.0

A longer horizon is a larger model¶

Repeat the day to make a longer horizon, and build the model without solving it. diagnostics gives its size:

def horizon(days: int) -> dict:
    hours = days * load.height
    repeated = pl.DataFrame({'snapshot': range(hours), 'value': days * load['value'].to_list()})
    return fleet | {'snapshot': range(hours), 'load': repeated, 'soc_initial': 0.0}


sizes = []
for days in (1, 3, 30, 365):
    with sps.build('tutorial/storage.yaml', horizon(days)) as model:
        size = model.diagnostics()
    sizes.append({'days': days, 'columns': size.columns, 'rows': size.rows})
print(pl.DataFrame(sizes))

shape: (4, 3)
┌──────┬─────────┬──────┐
│ days ┆ columns ┆ rows │
│ ---  ┆ ---     ┆ ---  │
│ i64  ┆ i64     ┆ i64  │
╞══════╪═════════╪══════╡
│ 1    ┆ 28      ┆ 24   │
│ 3    ┆ 84      ┆ 72   │
│ 30   ┆ 840     ┆ 720  │
│ 365  ┆ 10220   ┆ 8760 │
└──────┴─────────┴──────┘

Each day adds 28 columns and 24 rows. A real model has many more generators and hours than this one, and its solve gets slower as it grows.

Window by window¶

A rolling horizon solves the horizon in windows. Here three days stand in for a year, and the load is over hour:

days = pl.DataFrame({'hour': range(12), 'value': 3 * load['value'].to_list()})
hourly = fleet | {'load': days, 'soc_initial': 0.0}

EachWindow cuts hour into windows of three hours, and numbers the hours of each window from zero into snapshot. lookahead=0 means a window sees only the hours it keeps. carry copies the last soc of a window into the soc_initial of the next:

windows = sps.solve_over(
    'tutorial/storage.yaml',
    hourly,
    sps.EachWindow('hour', steps=3, lookahead=0, into='snapshot', carry={'soc_initial': 'soc'}),
)
print(windows.metrics.select('hour_start', 'columns', 'rows'))

shape: (4, 3)
┌────────────┬─────────┬──────┐
│ hour_start ┆ columns ┆ rows │
│ ---        ┆ ---     ┆ ---  │
│ i64        ┆ i64     ┆ i64  │
╞════════════╪═════════╪══════╡
│ 0          ┆ 21      ┆ 18   │
│ 3          ┆ 21      ┆ 18   │
│ 6          ┆ 21      ┆ 18   │
│ 9          ┆ 21      ┆ 18   │
└────────────┴─────────┴──────┘

Each window has 21 columns, however long the horizon is.

What a window cannot see¶

Solve the three days at once, and compare the cost of the two schedules. The objective of a window counts only that window, so add up what the windows kept:

def schedule_cost(sweep: sps.types.Sweep) -> float:
    output = sweep.primal('p').join(generators, on='generator')
    on = sweep.primal('on').join(generators, on='generator')
    return (output['value'] * output['cost']).sum() + (on['value'] * on['no_load_cost']).sum()


whole = sps.solve('tutorial/storage.yaml', horizon(3))
print(f'one solve {whole.objective:.1f}, windows {schedule_cost(windows):.1f}')

one solve 32550.0, windows 36461.1

The windows cost 12% more. To see why, look at the battery level at the end of each window:

print(windows.primal('soc').filter(pl.col('hour') % 3 == 2))

shape: (4, 2)
┌──────┬───────┐
│ hour ┆ value │
│ ---  ┆ ---   │
│ i64  ┆ f64   │
╞══════╪═══════╡
│ 2    ┆ 0.0   │
│ 5    ┆ 0.0   │
│ 8    ┆ 0.0   │
│ 11   ┆ 0.0   │
└──────┴───────┘

Every window ends with the battery empty. A window cannot see the hours after it, so it spends all the stored energy, and the next window needs gas sooner.

Lookahead¶

lookahead lets each window see more hours than it keeps, and solve them too. Try none, one and three:

rows = []
for lookahead in (0, 1, 3):
    sweep = sps.solve_over(
        'tutorial/storage.yaml',
        hourly,
        sps.EachWindow('hour', steps=3, lookahead=lookahead, into='snapshot', carry={'soc_initial': 'soc'}),
    )
    rows.append(
        {
            'lookahead': lookahead,
            'columns': sweep.metrics['columns'][0],
            'cost': round(schedule_cost(sweep), 1),
        }
    )
print(pl.DataFrame(rows))

shape: (3, 3)
┌───────────┬─────────┬─────────┐
│ lookahead ┆ columns ┆ cost    │
│ ---       ┆ ---     ┆ ---     │
│ i64       ┆ i64     ┆ f64     │
╞═══════════╪═════════╪═════════╡
│ 0         ┆ 21      ┆ 36461.1 │
│ 1         ┆ 28      ┆ 33255.6 │
│ 3         ┆ 42      ┆ 32550.0 │
└───────────┴─────────┴─────────┘

One hour of lookahead closes most of the gap, and three hours close all of it. Each hour also makes every window larger. The lookahead that is enough depends on the model: here, it is how far ahead the battery has to plan.

Next¶

These pages go further:


Sweeps and rolling horizons	every axis, `carry`, and how a sweep is read
Making a loop of solves faster	each window starts from the answer of the window before
Debugging a wrong answer	when a model solves and the number is wrong
Examples	every model in the repository