Change a model¶
The model is examples/dispatch.yaml — least-cost generation against a load
profile. Nothing here mutates it: every cell is a function of a spec and its
sources, so cells re-run in any order mean the same thing, and what you
carry out of the session is a file rather than a kernel.
Three loops, cheapest first:
- new numbers —
rebind, and the solver keeps the model it has loaded - more rows — the same math over a longer axis, which is still data
- new math — the spec is a
dict; patch it and re-run
and one verb that changes nothing: reading a built row back, for when the answer is wrong and the file looks right.
What none of them is: model.add_constraint(...). There is no Python API for
building a model here, on purpose — see the last cell.
import polars as pl
from IPython.display import Markdown
from math_spec import to_markdown, to_spec
import lpspec as lps
MODEL = '../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(MODEL))
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¶
$$\min \sum_{t \in \mathcal{T},\enspace 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\thinspace t \in \mathcal{T}$$
Variable domains¶
p
$$0 \le p_{t,g} \le \mathrm{p}^{\mathrm{max}}_{g} \qquad \forall\thinspace t \in \mathcal{T},\enspace g \in \mathcal{G} \thinspace:\thinspace \mathrm{p}^{\mathrm{max}}_{g} > 0$$
1. New numbers¶
build once, rebind per run of the cell. The declarations are untouched, so
the model HiGHS holds is untouched too: new costs go onto it in place, and the
matrix is never handed over a second time. The next solve begins from nothing
all the same — keep='solver', the default — because carrying the last
solve's work on costs the solver its presolve, which is a bet only you can
price (How much of the session a solve keeps). solve(keep='progress')
takes it.
bound = lps.build(MODEL, 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': bound.rebind({'cost': costs}).solve().objective})
sweep = pl.DataFrame(rows)
reused = bound.diagnostics()
print(f'{reused.loads} model loaded, {reused.solves} solves')
sweep
1 model loaded, 3 solves
| 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 model ever loaded.
That counter is the difference between "lpspec is slow" and "this loop is
rebuilding every time", and nothing about the answers depends on it.
The answers themselves are held to an equality: bound.rebind(x).solve() gives
what lps.solve(MODEL, sources | x) gives, always. Below, the last of the three
costs, solved from scratch. It is the oracle to reach for when a loop looks
wrong.
fresh = lps.solve(MODEL, sources | {'cost': costs}).objective
rebound = sweep.filter(pl.col('gas_cost') == 90.0).item(0, 'objective')
print(f'rebound {rebound:,.1f} — fresh build {fresh:,.1f}')
rebound 9,900.0 — fresh build 9,900.0
2. More rows¶
A longer horizon is not a different model — it is a longer table plus the
index to match. The same move grows a Benders cut family one cut at a time
(examples/benders/run.py), which is what makes "add a constraint" a data
question far more often than it looks.
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 = bound.rebind({'snapshot': index, 'load': horizon}).solve().primal('p')
grown = bound.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
| 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 that model was
loaded from scratch and solved cold — the fast path is what changes, never the
answer. Read the counter, not the clock.
p comes back tidy and in label order, so nothing here has to sort:
schedule.pivot(on='generator', index='snapshot', values='value') is the wide
view if you would rather read a schedule than a table of rows.
3. New math¶
to_dict() is the model as data, and every verb takes a dict as readily as a
path. So an edit is a key, and the round trip through to_spec re-validates
it — the language's load-time errors are this notebook's error messages.
Below: a ramp limit on gas, which needs a parameter as well as a constraint.
spec = to_spec(MODEL).to_dict()
spec['parameters']['ramp_max'] = {'dims': ['generator']}
spec['constraints']['ramp_up'] = {
'foreach': ['snapshot', 'generator'],
'expression': 'p - shift(p, over=snapshot, offset=1) <= ramp_max',
}
ramp_max = pl.DataFrame({'generator': GENERATORS, 'value': [100.0, 100.0, 20.0]})
base = lps.solve(MODEL, sources).objective
ramped = lps.solve(spec, sources | {'ramp_max': ramp_max}).objective
pl.DataFrame({'model': ['dispatch', 'dispatch + ramp limit'], 'objective': [base, ramped]})
| model | objective |
|---|---|
| str | f64 |
| "dispatch" | 6600.0 |
| "dispatch + ramp limit" | 8400.0 |
The limit binds: gas cannot reach the evening peak in one step, so it starts climbing early and free wind is curtailed to make room for it.
The math re-renders from the patched spec, which is the check that the edit says what you meant:
Markdown(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},\enspace 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\thinspace t \in \mathcal{T}$$
ramp_up
$$p_{t,g} - p_{t - 1,g} \le \mathrm{ramp\_max}_{g} \qquad \forall\thinspace t \in \mathcal{T},\enspace g \in \mathcal{G}$$
Variable domains¶
p
$$0 \le p_{t,g} \le \mathrm{p}^{\mathrm{max}}_{g} \qquad \forall\thinspace t \in \mathcal{T},\enspace g \in \mathcal{G} \thinspace:\thinspace \mathrm{p}^{\mathrm{max}}_{g} > 0$$
An edit the language refuses is refused before any data is bound — check
parses, resolves and lowers, and nothing else runs:
typo = {
**spec,
'constraints': {
**spec['constraints'],
'peak': {'foreach': ['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. Dimensions: ['generator', 'snapshot'] 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¶
An edit the language refuses never reaches the data. An edit it accepts can
still build something other than what the file appears to say, because absence
is how this language masks: 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.
Gas is retired below the cheapest way there is, one number, and the model stops having an answer.
retired = sources | {'p_max': pl.DataFrame({'generator': GENERATORS, 'value': [80.0, 40.0, 0.0]})}
short = lps.build(MODEL, 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 against a load three generators cover comfortably, and the file
still reads sum(p, over=generator) == load over all three. row gives the
row the build produced at one coordinate, off the built model and with no
solve — so it answers on a model that never reached a solver:
fleet = lps.build(MODEL, 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 not zero in the second row, it is missing from it: six columns
went with the where, and power_balance is left asking two generators for a
load of 180. The line is linopy's shape for
the same job, with the row's own coordinate added to it.
That is the class of fault no solver output names, because the solver was
handed a model perfectly consistent with itself. Where a mask takes every row
of a declaration, row says so at the coordinate and diagnostics().omissions
counts them — the only report a constraint that built nothing gets.
What leaves the session¶
The spec, as the file you diff against examples/dispatch.yaml and commit —
not this notebook, and not a pickle of the kernel. A model built as a dict still gets
a file, which is the whole point of the round trip:
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:
foreach:
- 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:
foreach:
- snapshot
expression: sum(p, over=generator) == load
ramp_up:
foreach:
- snapshot
- generator
expression: p - shift(p, over=snapshot, offset=1) <= ramp_max
objective:
sense: minimize
expression: sum(p * cost)
description: total cost of generation over the horizon
Where this goes next¶
The verbs a linopy reader reaches for — fix, relax, "remove that constraint"
— are these same three loops aimed at a different question, and they have their
own page: Fix, relax, remove, which also names what has no
spelling here at all.
What the loops buy is the other side of that trade: a cell that cannot half-apply, a session whose output is a diff, and re-run order that cannot change what the model means.