Skip to content

Change a model

The dispatch model from Run a model, changed three ways, cheapest first. Every block is a function of a spec (the YAML, as a dict) and its sources, so blocks re-run in any order mean the same thing.

  1. New numbers: update, and the solver keeps the model it has loaded.
  2. More rows: the same math over a longer axis.
  3. New math: patch the dict and re-run.

Every block on this page runs when the site is built, and what you see under it is what it printed on this commit. A block that raises fails the build.

import polars as pl
from mathspec import to_markdown, to_spec

import specsolve as sps

SPEC = 'examples/dispatch.yaml'
GENERATORS = ['wind', 'solar', 'gas']

sources = {
    'snapshot': pl.DataFrame({'snapshot': range(6)}),
    'generator': pl.DataFrame({'generator': GENERATORS}),
    'p_max': pl.DataFrame({'generator': GENERATORS, 'value': [80.0, 40.0, 200.0]}),
    'cost': pl.DataFrame({'generator': GENERATORS, 'value': [0.0, 0.0, 60.0]}),
    'load': pl.DataFrame({'snapshot': range(6), 'value': [90.0, 120.0, 150.0, 180.0, 140.0, 100.0]}),
}

print(to_markdown(SPEC))

Least-cost dispatch of a generator fleet against an hourly load.

Sets

Symbol Meaning
\(\mathcal{T}\) index \(t\) — snapshot — dispatch periods
\(\mathcal{G}\) index \(g\) — generator — generating units

Parameters

Symbol Meaning
\(\mathrm{p}^{\mathrm{max}}\) p_max over \(\mathcal{G}\) — installed capacity
\(\mathrm{load}\) load over \(\mathcal{T}\) — demand to be met
\(\mathrm{cost}\) cost over \(\mathcal{G}\) — marginal cost

Variables

Symbol Meaning
\(p\) p over \(\mathcal{T} \times \mathcal{G}\) — output of a generator in a snapshot

Upright is what the data supplies — a parameter such as \(\mathrm{p}^{\mathrm{max}}\), a coordinate map, a label — and italic is what the solver chooses, such as \(p\). An index is italic too, being what a quantifier chooses, and a set is script.

Objective

\[ \min \sum_{t \in \mathcal{T},\ g \in \mathcal{G}} p_{t,g} \cdot \mathrm{cost}_{g} \]

Subject to

power_balance

\[ \sum_{g \in \mathcal{G}} p_{t,g} = \mathrm{load}_{t} \qquad \forall\, t \in \mathcal{T} \]

Variable domains

p

\[ 0 \le p_{t,g} \le \mathrm{p}^{\mathrm{max}}_{g} \qquad \forall\, t \in \mathcal{T},\ g \in \mathcal{G} \,:\, \mathrm{p}^{\mathrm{max}}_{g} > 0 \]

1. New numbers

Build once, then update per run. New costs go onto the model HiGHS holds, and the matrix is never handed over twice.

model = sps.build(SPEC, sources)

rows = []
for gas_cost in (40.0, 60.0, 90.0):
    costs = pl.DataFrame({'generator': GENERATORS, 'value': [0.0, 0.0, gas_cost]})
    rows.append({'gas_cost': gas_cost, 'objective': model.update({'cost': costs}).solve().objective})

sweep = pl.DataFrame(rows)
reused = model.diagnostics()

print(f'{reused.loads} model loaded, {reused.solves} solves')
print(sweep)

1 model loaded, 3 solves
shape: (3, 2)
┌──────────┬───────────┐
│ gas_cost ┆ objective │
│ ---      ┆ ---       │
│ f64      ┆ f64       │
╞══════════╪═══════════╡
│ 40.0     ┆ 4400.0    │
│ 60.0     ┆ 6600.0    │
│ 90.0     ┆ 9900.0    │
└──────────┴───────────┘

loads is 1 against solves of 3: three answers, one load. model.update(x).solve() gives what sps.solve(SPEC, sources | x) gives, always. The next block solves the last cost from scratch to show it.

fresh = sps.solve(SPEC, sources | {'cost': costs}).objective
updated = sweep.filter(pl.col('gas_cost') == 90.0).item(0, 'objective')

print(f'updated {updated:,.1f} — fresh build {fresh:,.1f}')

updated 9,900.0 — fresh build 9,900.0

2. More rows

A longer horizon is a longer table plus the index to match.

horizon = pl.DataFrame(
    {
        'snapshot': range(12),
        'value': [90.0, 120.0, 150.0, 180.0, 140.0, 100.0, 95.0, 130.0, 160.0, 190.0, 150.0, 110.0],
    }
)

index = pl.DataFrame({'snapshot': range(12)})
schedule = model.update({'snapshot': index, 'load': horizon}).solve().primal('p')
grown = model.diagnostics()

print(f'{schedule.height} rows of p now, and {grown.loads} loads over {grown.solves} solves')
print(schedule.head())

36 rows of p now, and 2 loads over 4 solves
shape: (5, 3)
┌──────────┬───────────┬───────┐
│ snapshot ┆ generator ┆ value │
│ ---      ┆ ---       ┆ ---   │
│ i64      ┆ str       ┆ f64   │
╞══════════╪═══════════╪═══════╡
│ 0        ┆ wind      ┆ 50.0  │
│ 0        ┆ solar     ┆ 40.0  │
│ 0        ┆ gas       ┆ 0.0   │
│ 1        ┆ wind      ┆ 80.0  │
│ 1        ┆ solar     ┆ 40.0  │
└──────────┴───────────┴───────┘

loads is 2 now: new coordinates renumber the columns, so this model was loaded from scratch. The answer is the same either way. schedule.pivot(on='generator', index='snapshot', values='value') is the wide view.

3. New math

to_dict() is the spec as data, and every verb takes a dict. An edit is a key, and to_spec validates it again. Below, a ramp limit on gas, which needs a parameter as well as a constraint.

spec = to_spec(SPEC).to_dict()
spec['parameters']['ramp_max'] = {'dims': ['generator']}
spec['constraints']['ramp_up'] = {
    'dims': ['snapshot', 'generator'],
    'expression': 'p - shift(p, along=snapshot, offset=1) <= ramp_max',
}

ramp_max = pl.DataFrame({'generator': GENERATORS, 'value': [100.0, 100.0, 20.0]})
base = sps.solve(SPEC, sources).objective
ramped = sps.solve(spec, sources | {'ramp_max': ramp_max}).objective

print(pl.DataFrame({'model': ['dispatch', 'dispatch + ramp limit'], 'objective': [base, ramped]}))

shape: (2, 2)
┌───────────────────────┬───────────┐
│ model                 ┆ objective │
│ ---                   ┆ ---       │
│ str                   ┆ f64       │
╞═══════════════════════╪═══════════╡
│ dispatch              ┆ 6600.0    │
│ dispatch + ramp limit ┆ 8400.0    │
└───────────────────────┴───────────┘

The limit binds: gas starts climbing early, and free wind is curtailed to make room. The math re-renders from the patched spec:

print(to_markdown(spec, legend=False, numbered=False))

Least-cost dispatch of a generator fleet against an hourly load.

Objective

\[ \min \sum_{t \in \mathcal{T},\ g \in \mathcal{G}} p_{t,g} \cdot \mathrm{cost}_{g} \]

Subject to

power_balance

\[ \sum_{g \in \mathcal{G}} p_{t,g} = \mathrm{load}_{t} \qquad \forall\, t \in \mathcal{T} \]

ramp_up

\[ p_{t,g} - p_{t - 1,g} \le \mathrm{ramp\_max}_{g} \qquad \forall\, t \in \mathcal{T},\ g \in \mathcal{G} \]

Variable domains

p

\[ 0 \le p_{t,g} \le \mathrm{p}^{\mathrm{max}}_{g} \qquad \forall\, t \in \mathcal{T},\ g \in \mathcal{G} \,:\, \mathrm{p}^{\mathrm{max}}_{g} > 0 \]

An edit the language refuses is refused before any data is attached:

typo = {
    **spec,
    'constraints': {
        **spec['constraints'],
        'peak': {'dims': ['snapshot'], 'expression': 'sum(p, over=generators) <= load'},
    },
}

try:
    sps.check(typo)
except sps.errors.LanguageError as exc:
    print(exc)

Constraint 'peak': sum(over=generators) does not name a declared dimension. Did you mean 'generator'?
Declare 'generators' under 'dimensions:', or fix the typo — an unknown dimension makes sum() a silent no-op rather than an error.

What leaves the session

The spec, as a file you can diff and commit:

print(to_spec(spec).to_yaml())

version: 0
description: Least-cost dispatch of a generator fleet against an hourly load.
dimensions:
  snapshot:
    dtype: int
    description: dispatch periods
  generator:
    dtype: str
    description: generating units
parameters:
  p_max:
    dims:
    - generator
    dtype: float
    description: installed capacity
  load:
    dims:
    - snapshot
    dtype: float
    description: demand to be met
  cost:
    dims:
    - generator
    dtype: float
    description: marginal cost
  ramp_max:
    dims:
    - generator
    dtype: float
variables:
  p:
    dims:
    - snapshot
    - generator
    where: p_max > 0
    bounds:
      lower: 0.0
      upper: p_max
    domain: continuous
    absence: undefined
    description: output of a generator in a snapshot
constraints:
  power_balance:
    dims:
    - snapshot
    expression: sum(p, over=generator) == load
  ramp_up:
    dims:
    - snapshot
    - generator
    expression: p - shift(p, along=snapshot, offset=1) <= ramp_max
objective:
  sense: minimize
  expression: sum(p * cost)
  description: total cost of generation over the horizon

Where next
Sweep a model the next tutorial: one model once per scenario, then window by window
Warm-starting a re-solve keep the solver's work between solves, not only the model
Fixing, relaxing and removing linopy's fix, relax and remove_constraints, as these three loops
Debugging a wrong answer read the row a build produced, when the file looks right and the answer is not