Change a model¶

The dispatch model from Run a model, changed three ways, cheapest first. Every cell is a function of a spec (the YAML, as a dict) and its sources, so cells 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.

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

In [1]:

import polars as pl
from IPython.display import Markdown
from math_spec import to_markdown, to_spec

import lpspec as lps

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]}),
}

Markdown(to_markdown(SPEC))

import polars as pl from IPython.display import Markdown from math_spec import to_markdown, to_spec import lpspec as lps 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]}), } Markdown(to_markdown(SPEC))

Out[1]:

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 model is given — 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¶

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

Subject to¶

power_balance

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

Variable domains¶

p

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

In [2]:

model = lps.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')
sweep

model = lps.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') sweep

1 model loaded, 3 solves

Out[2]:

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 lps.solve(SPEC, sources | x) gives, always. The next cell solves the last cost from scratch to show it.

In [3]:

fresh = lps.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}')

fresh = lps.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.

In [4]:

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')
schedule.head()

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') schedule.head()

36 rows of p now, and 2 loads over 4 solves

Out[4]:

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.

In [5]:

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 = lps.solve(SPEC, sources).objective
ramped = lps.solve(spec, sources | {'ramp_max': ramp_max}).objective

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

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 = lps.solve(SPEC, sources).objective ramped = lps.solve(spec, sources | {'ramp_max': ramp_max}).objective pl.DataFrame({'model': ['dispatch', 'dispatch + ramp limit'], 'objective': [base, ramped]})

Out[5]:

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:

In [6]:

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

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

Out[6]:

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

Objective¶

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

Subject to¶

power_balance

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

ramp_up

math
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

math
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:

In [7]:

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

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

typo = { **spec, 'constraints': { **spec['constraints'], 'peak': {'dims': ['snapshot'], 'expression': 'sum(p, over=generators) <= load'}, }, } try: lps.check(typo) except lps.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:

In [8]:

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

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

retired = sources | {'p_max': pl.DataFrame({'generator': GENERATORS, 'value': [80.0, 40.0, 0.0]})} short = lps.build(SPEC, retired) answer = short.solve() print(f'{answer.status} / {answer.termination_condition}') try: answer.primal('p') except lps.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:

In [9]:

fleet = lps.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}')

fleet = lps.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:

In [10]:

print(to_spec(spec).to_yaml())

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. How much of the session a solve keeps is the one choice made for you here: keep='solver'.