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.

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

A fourth section reads a built row back, for when the answer is wrong and the file looks right.

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.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.

4. When the answer is wrong¶

p is declared where: "p_max > 0", so a capacity of zero does not park a generator at zero. It deletes the column and every term that referenced it. Below, gas is retired with one number:

retired = sources | {'p_max': pl.DataFrame({'generator': GENERATORS, 'value': [80.0, 40.0, 0.0]})}
short = sps.build(SPEC, retired)
answer = short.solve()

print(f'{answer.status} / {answer.termination_condition}')
try:
    answer.primal('p')
except sps.NoSolutionError as exc:
    print(exc)

warning / infeasible
cannot read the primal of 'p': the solve terminated 'infeasible' (Infeasible), so there are no values to read. Test `has_primal` first. This raises rather than returning, because the solver hands back a full-length vector of zeros either way and it is indistinguishable from an answer.

Infeasible, while the file still reads sum(p, over=generator) == load over all three generators. row gives the row the build produced, with no solve:

fleet = sps.build(SPEC, sources)

print(fleet.row('power_balance', snapshot=3))
print(short.row('power_balance', snapshot=3))
print(f'{fleet.diagnostics().columns} columns became {short.diagnostics().columns}')

power_balance[snapshot=3]: +1 p[3, wind] +1 p[3, solar] +1 p[3, gas] == 180
power_balance[snapshot=3]: +1 p[3, wind] +1 p[3, solar] == 180
18 columns became 12

p[3, gas] is missing from the second row: six columns went with the where, and power_balance asks two generators for a load of 180. The line is linopy's shape. No solver output names this fault. Where a mask takes every row of a declaration, diagnostics().omissions counts them.

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¶

Fix, relax, remove spells fix, relax and "remove that constraint" as these three loops. keep= is the one choice made for you here: 'solver'.