Energiesysteemmodellering met lineair programmeren

Energiesysteemmodellen proberen een praktische vraag precies genoeg te maken om ermee te kunnen rekenen. Niet: “hoe ziet een goed energiesysteem eruit?”, maar bijvoorbeeld:

Welke mix van wind en gas kan in elk uur genoeg elektriciteit leveren, tegen de laagste kosten?

Dat klinkt meteen als een optimalisatieprobleem. Je hebt getallen die vaststaan, keuzes die je mag maken, eisen waaraan het systeem moet voldoen, en een doel dat je zo goed mogelijk wilt maken. In veel energiemodellen worden die onderdelen vertaald naar een lineair programma:

Parameters: getallen die de gebruiker vastlegt, zoals vraag, windfactoren en kosten.
Decision variables, of beslissingsvariabelen: de onbekenden die de solver mag kiezen.
Constraints, of beperkingen: regels waaraan elke oplossing moet voldoen.
Objective, of doelfunctie: wat je zo klein of groot mogelijk wilt maken.

Een handige manier om dit uit elkaar te houden is:

Type	Wie bepaalt dit?
Parameters	de gebruiker
Decision variables	de solver
Constraints	het model
Objective	het model

Het woord “lineair” betekent hier dat de variabelen alleen in rechte combinaties voorkomen. Zowel de constraints als de objective bestaan uit gewogen optellingen van variabelen. Elke variabele wordt hoogstens vermenigvuldigd met een constante. Daardoor blijven de grenzen in de figuur rechte lijnen.

Geen $W^2$ , geen $1/G$ , geen “als-dan” sprongen dus. Die beperking lijkt streng, maar ze is precies waarom zulke modellen zo krachtig zijn: lineaire programma’s zijn goed oplosbaar, ook wanneer ze duizenden variabelen en constraints hebben.

Een mini-energiesysteem

We bouwen eerst een heel klein model. Stel dat een regio elektriciteit kan leveren met twee technologieën:

windvermogen, aangeduid met $W$ ;
gasvermogen, aangeduid met $G$ .

$W$ en $G$ zijn de eerste twee decision variables van ons model. Beide variabelen stellen capaciteit voor. Als het model kiest voor $W = 80$ , dan betekent dat: bouw 80 GW windcapaciteit. Als het kiest voor $G = 40$ , dan betekent dat: houd 40 GW gasvermogen beschikbaar.

De eerste regels zijn logisch:

W \ge 0 \qquad G \ge 0

Negatieve windmolens of negatieve gascentrales bestaan niet. Het model mag dus alleen punten kiezen in het eerste kwadrant: rechts en boven nul.

Nu voegen we twee representatieve uren toe. In het eerste uur waait het weinig. Wind kan dan maximaal 25 procent van zijn geïnstalleerde capaciteit leveren. De vraag is 60 GW:

0.25W + G \ge 60

Merk op dat we hier bewust vereenvoudigen. We modelleren nog geen aparte productie per uur, maar alleen of er genoeg beschikbare capaciteit is om de vraag te kunnen dekken. In een echt energiemodel zijn capaciteit en productie verschillende variabelen. Voor deze eerste kennismaking houden we ze samen, zodat de geometrie van lineair programmeren zichtbaar blijft.

In het tweede uur waait het beter. Wind kan dan maximaal 60 procent van zijn capaciteit leveren, maar de vraag is ook hoger: 90 GW.

0.60W + G \ge 90

Dit is al een echt energiemodel in miniatuur. Wind is in dit voorbeeld per GW goedkoper, maar is niet altijd volledig beschikbaar. Gas is flexibel, maar duurder. Het model moet een combinatie vinden die in beide situaties genoeg elektriciteit kan leveren.

De getallen $60$ , $90$ , $0.25$ en $0.60$ zijn parameters. Wij kiezen ze als invoer voor het model. De solver mag ze niet veranderen; hij mag alleen waarden zoeken voor de decision variables $W$ en $G$ .

De doelfunctie

Een optimalisatiemodel heeft niet alleen regels nodig, maar ook een richting. Hier willen we kosten minimaliseren. Stel dat een extra GW wind jaarlijks 20 kost en een extra GW gas jaarlijks 45. Dan wordt de doelfunctie:

\min \quad 20W + 45G

Alles bij elkaar krijgen we:

\begin{aligned} \min \quad &20W + 45G \\ \text{zodat} \quad &0.25W + G \ge 60 \\ &0.60W + G \ge 90 \\ &W \ge 0,\quad G \ge 0 \end{aligned}

Dat is de hele puzzel. Het mooie is dat je hem kunt tekenen.

Op dit moment hebben we nog geen windmolens gebouwd of elektriciteit geproduceerd. We hebben alleen een verzameling toegestane punten in een tweedimensionale ruimte beschreven. Elk punt stelt één mogelijke keuze voor wind- en gascapaciteit voor. De taak van de solver is nu verrassend eenvoudig: zoek binnen die ruimte het punt met de laagste kosten.

In dit voorbeeld zijn de decision variables alleen $W$ en $G$ . In een echt energiemodel zijn dat er vaak honderdduizenden: dispatch per uur, batterij-inhoud, lijnstromen, nieuwe capaciteit, en meer. De solver kiest al die decision variables tegelijk.

Wat je in de figuur ziet

De visualisatie is opgebouwd in vier stappen. Dat is bewust: bij lineair programmeren wil je eerst zien welke oplossingen toegestaan zijn, en pas daarna welke toegestane oplossing het goedkoopst is.

1. Constraints

Elke constraint is een lijn die een deel van het vlak uitsluit. De lijn

0.25W + G = 60

is precies de grens tussen “genoeg beschikbare capaciteit in het windarme uur” en “te weinig beschikbare capaciteit in het windarme uur”. Omdat de constraint $\ge 60$ is, is alles onder die lijn uitgesloten. Hetzelfde geldt voor het windrijke uur, maar met een andere windfactor en een andere vraag. De twee kleuren laten dus twee aparte redenen zien waarom een punt niet mag.

2. Overlap

Wat na alle uitsluitingen overblijft, is de feasible region: alle combinaties van $W$ en $G$ die aan beide constraints tegelijk voldoen. Elk punt in dat overblijvende gebied is technisch toegestaan. Het model mag dus bijvoorbeeld veel te veel capaciteit bouwen. Dat is niet verboden door de constraints. Het is alleen duur.

De sliders voor vraag en windfactoren veranderen deze ruimte direct. Hogere vraag duwt een constraint omhoog. Een lagere windfactor maakt wind minder effectief, waardoor de lijn vlakker wordt en er meer wind- of gascapaciteit nodig is.

De feasible region is in dit voorbeeld eigenlijk onbegrensd: je kunt altijd meer wind en gas bouwen dan nodig. De figuur knipt die ruimte af op een zichtbaar venster, zodat je de interessante randen en hoekpunten kunt zien.

Probeer in de visualisatie eens de windfactor van het windarme uur naar 0.1 te schuiven. De constraint wordt vlakker en de feasible region verschuift zichtbaar. Zonder extra algebra zie je meteen dat het systeem meer capaciteit nodig heeft.

3. Hoekpunten

In de derde stap markeren we de kandidaat-hoekpunten van de feasible region. Die punten zijn niet vooraf gekozen. Ze worden berekend als snijpunten van de twee constraintgrenzen en de assen $W = 0$ en $G = 0$ . Die assen tellen mee omdat ook $W \ge 0$ en $G \ge 0$ constraints zijn.

Dit is de belangrijkste intuitiesprong achter lineair programmeren:

Als een lineair programma een optimum heeft, dan is er altijd minstens één optimum op een hoekpunt van het toegestane gebied.

Soms is niet één punt uniek het beste, maar is een hele rand even goed. Ook dan ligt er minstens één optimale oplossing op een hoekpunt. Daarom hoeft een solver niet elk punt in het groene vlak te proberen. Hij hoeft slim door de hoekpunten te bewegen. Het simplex-algoritme doet precies dat op grote schaal. In deze blog gebruiken we niet het volledige tableau-algoritme; we bouwen eerst de hoekpunt-logica zelf na.

Andere moderne algoritmen, zoals interior-point methoden, volgen een andere route door de oplossingsruimte. Toch gebruiken ze dezelfde lineaire beschrijving van het probleem en vinden ze uiteindelijk hetzelfde optimum.

4. Kosten

De gestippelde lijn is een gelijke-kostenlijn. Op zo’n lijn hebben alle punten dezelfde waarde van

20W + 45G

De kosten-sliders verschijnen pas in deze stap, omdat kosten de feasible region niet veranderen. Ze veranderen alleen de helling van de objective-lijn en dus welk hoekpunt het goedkoopst is. Als wind goedkoper wordt ten opzichte van gas, draait de gelijke-kostenlijn zichtbaar mee. Als een ander hoekpunt goedkoper wordt, verspringt het optimum.

Standaard raakt de kostlijn het optimale hoekpunt. Met “kostlijn schuiven” kun je hem alleen demonstratief verplaatsen: lagere kosten liggen buiten de feasible region, hogere kosten zijn toegestaan maar niet optimaal.

Een goede test voor je intuïtie: zet de gaskosten lager dan de windkosten. Kun je voorspellen welk hoekpunt optimaal wordt voordat je naar de tabel kijkt?

De hoekpunten zelf vinden

In twee dimensies kun je de kern van het idee met gewone Python namaken. Elke hoek ontstaat waar twee grenslijnen elkaar snijden. Daarna filter je de punten die niet feasible zijn, bereken je de kosten, en kies je het goedkoopste punt.

from itertools import combinations

constraints = [
    # a*W + b*G >= rhs
    (0.25, 1.0, 60.0, "windarm uur"),
    (0.60, 1.0, 90.0, "windrijk uur"),
    (1.00, 0.0, 0.0, "W >= 0"),
    (0.00, 1.0, 0.0, "G >= 0"),
]

def intersect(c1, c2):
    a1, b1, rhs1, _ = c1
    a2, b2, rhs2, _ = c2
    det = a1 * b2 - a2 * b1
    if abs(det) < 1e-9:
        return None
    W = (rhs1 * b2 - rhs2 * b1) / det
    G = (a1 * rhs2 - a2 * rhs1) / det
    return W, G

def feasible(point):
    W, G = point
    return all(a * W + b * G >= rhs - 1e-9 for a, b, rhs, _ in constraints)

def cost(point):
    W, G = point
    return 20 * W + 45 * G

points = []

for c1, c2 in combinations(constraints, 2):
    point = intersect(c1, c2)
    if point is not None and feasible(point):
        points.append(point)

best = min(points, key=cost)

print(best)
print(cost(best))

De uitkomst is ongeveer:

(85.71428571428571, 38.57142857142857)
3450.0

Dus het goedkope systeem in dit mini-model bouwt ongeveer 85.7 GW wind en 38.6 GW gasvermogen. Niet omdat gas “slecht” is of wind “goed”, maar omdat precies deze combinatie op beide representatieve uren genoeg levert tegen de laagste kosten.

Tot nu toe hebben wij steeds over wind en gas gesproken. Voor de solver bestaan die woorden niet. De solver ziet geen windmolens, geen gascentrales en geen elektriciteit. Hij ziet alleen variabelen, lineaire constraints en een objective. De betekenis van die variabelen komt volledig uit het model dat wij hebben opgesteld.

Waarom dit al op PyPSA lijkt

PyPSA doet hetzelfde soort werk, maar met veel meer structuur. PyPSA is niet zelf de solver. Het is een modelgenerator: je beschrijft een energiesysteem, PyPSA zet dat om in een lineair optimalisatieprobleem, en een externe solver zoals HiGHS lost dat probleem vervolgens op.

In plaats van zelf lijsten met constraints te schrijven, beschrijf je een netwerk:

een bus is een knooppunt waar elektriciteit samenkomt;
een generator is een technologie die aan zo’n bus kan leveren;
een snapshot is een tijdstap, bijvoorbeeld een uur;
per snapshot moet de energiebalans kloppen;
de objective telt investeringskosten, operationele kosten, of allebei.

Ons mini-model heeft eigenlijk al PyPSA-vormen:

In het mini-model	In PyPSA-taal
$W$	uitbreidbare windcapaciteit
$G$	uitbreidbare gascapaciteit
windfactoren 0.25 en 0.60	tijdsafhankelijke beschikbaarheid
vraag 60 en 90	load per snapshot
$0.25W + G \ge 60$	capaciteits-/voorzieningsconstraint
$20W + 45G$	objective met capaciteitskosten

Strikt genomen is dit mini-model nog eenvoudiger dan PyPSA. We laten productievariabelen weg en doen alsof beschikbare capaciteit direct voldoende is om de vraag te dekken. PyPSA maakt die stap expliciet. Het kiest niet alleen hoeveel wind- en gascapaciteit er gebouwd wordt, maar ook hoeveel elke technologie produceert in elke tijdstap. Voor het windarme uur zou PyPSA onder de motorkap eerder iets maken als:

\begin{aligned} p_{\text{wind},1} + p_{\text{gas},1} &= 60 \\ p_{\text{wind},1} &\le 0.25W \\ p_{\text{gas},1} &\le G \end{aligned}

Dat verschil wordt belangrijk zodra je marginale kosten, brandstofkosten of dispatch wilt uitleggen. In ons mini-model staan alleen capaciteitskosten in de objective:

20W + 45G

Met dispatchkosten zou daar bijvoorbeeld iets bijkomen als:

20W + 45G + 70p_{\text{gas},1} + 70p_{\text{gas},2}

Voor de eerste intuïtie is het compacte capaciteitsmodel genoeg. Voor echte PyPSA-modellen is het belangrijk om te onthouden dat capaciteit en productie twee verschillende soorten variabelen zijn.

Een schets in PyPSA-achtige code zou er conceptueel zo uitzien:

import pypsa

n = pypsa.Network()
n.set_snapshots(["windarm", "windrijk"])

n.add("Bus", "elektriciteit")

n.add(
    "Load",
    "vraag",
    bus="elektriciteit",
    p_set=[60, 90],
)

n.add(
    "Generator",
    "wind",
    bus="elektriciteit",
    p_nom_extendable=True,
    capital_cost=20,
    p_max_pu=[0.25, 0.60],
)

n.add(
    "Generator",
    "gas",
    bus="elektriciteit",
    p_nom_extendable=True,
    capital_cost=45,
)

n.optimize()

De namen zijn anders, maar de logica is dezelfde. PyPSA maakt onder de motorkap de variabelen, constraints en objective aan. Daarna geeft het model die lineaire puzzel door aan een solver.

Alles wat we hierboven met de hand hebben gedaan, doet PyPSA dus automatisch:

Met de hand	In PyPSA
$W$ en $G$ kiezen	decision variables maken
constraints opschrijven	constraints genereren
`cost(point)` berekenen	objective opbouwen
`min(points, key=cost)`	solver aanroepen

Waar de echte complexiteit begint

Het mini-model is bewust te klein. Echte energiesysteemmodellen voegen veel lagen toe:

meerdere regio’s en transmissielijnen;
veel meer uren, vaak een heel jaar;
opslag met laden, ontladen en energie-inhoud;
CO2-limieten;
bestaande capaciteit naast nieuwe investeringen;
brandstofkosten, rendementen en operationele beperkingen.

Maar de taal blijft herkenbaar. Je blijft vragen:

Welke variabelen mag het model kiezen?
Welke constraints mogen nooit geschonden worden?
Welke objective geeft richting aan de keuze?

Als je die drie vragen kunt beantwoorden, kun je een energiemodel lezen. En als je de feasible region hierboven begrijpt, heb je de kern van lineair programmeren al in handen: een model zoekt niet magisch naar “het beste energiesysteem”, maar beweegt door een geometrische ruimte van mogelijke oplossingen totdat het een goedkoop toegestaan hoekpunt vindt. Hoe een solver dat vervolgens efficiënt doet, is weer een interessant onderwerp op zichzelf.