Companion code for:
Till Holmes, Sverre Stefanussen Foslie, Stian Backe and Hanne Kauko. Thermal Energy Storage and District Heating as Flexibility Enablers for Industrial Heat Decarbonization.
A linear investment-and-dispatch optimisation model, built on EnergyModelsX, that couples industrial process heat with a district heating network and the electricity grid for a region in Germany. It minimises total system cost across five five-year investment periods under a tightening CO2 cap, and is used to quantify how thermal energy storage (TES) and industry–district-heating heat exchange change the cost-optimal decarbonization pathway.
Three industry sectors are represented, each with low-, medium- and high-temperature heat demand, gas and electric boilers, heat pumps and thermal storage. Each investment period is one full year at hourly resolution (8760 operational periods).
- Julia 1.12 (download). The committed
Manifest.tomlwas resolved with 1.12.6; that is what reproduces the environment exactly. On an older Julia, deleteManifest.tomland re-resolve. - ~16 GB RAM. Even the smallest configuration builds an LP with roughly 4.5 million rows and uses several GB.
- A solver. HiGHS is open source and installed automatically. Gurobi is optional and commercial — see Solvers.
julia --project -e 'using Pkg; Pkg.instantiate()'Then run the smoke-test configuration, which takes about 40 seconds and exists purely to confirm your setup works:
julia --project run_case.jl config.quick.tomlResults land in results/quick/. Those numbers are not scientific — the quick config
shortens the year to 30 days. Do not start with config.toml; read
Solvers first.
The published results were produced with Gurobi. This matters more than it usually does, because the model is large: one industry sector over a single full-year investment period already yields an LP with ~4.5 million rows and ~3.4 million columns.
Measured on an Apple M-series laptop:
| Configuration | Sectors × periods × hours | LP size | HiGHS | Gurobi |
|---|---|---|---|---|
config.quick.toml |
1 × 1 × 720 | 273 k rows | ~40 s | not measured |
| single sector, full year | 1 × 1 × 8760 | 4.5 M rows | did not converge in 39 min | not measured |
config.toml (paper) |
3 × 5 × 8760 | ~15× the above | impractical | used for the paper |
Note the middle row — that is the smallest scientifically meaningful configuration, and
it did not finish. Shortening the year from 8760 to 720 hours is the difference between
40 seconds and "come back tomorrow", which is why config.quick.toml truncates the
horizon: without it there is no configuration you can use to check your setup quickly.
With HiGHS:
- Run the four cases individually (
run_case.jl), not viarun_all_cases.jl. A single failure or restart then does not cost you every case. - Expect long solves.
run_case.jlprints solver progress so you can see it is working.
Gurobi is deliberately not a dependency of this project, so that readers without a licence get a small, fast install. If you have a licence:
julia --project -e 'using Pkg; Pkg.add("Gurobi")'Then set name = "Gurobi" under [solver] in your config. The code loads it on demand
and gives an actionable error if it is missing.
The paper reports four scenarios, crossing TES availability with the industry–district heating connection (Table 2):
| Case | tes_enabled |
dh_connection |
|---|---|---|
| Base (no TES, separated) | false |
false |
| TES, separated | true |
false |
| No TES, connected | false |
true |
| TES, connected | true |
true |
Disabling TES sets its investment cost prohibitively high rather than removing the nodes, so the variable structure stays identical across cases. Disabling the connection caps the district-heating connection capacity at zero.
Run one case by editing config.toml, or all four with:
julia --project run_all_cases.jlEverything about a run lives in a TOML file, so a run is fully described by one readable
artifact you can commit next to its results. See config.toml for the
paper setup and config.quick.toml for the smoke test. Keys:
| Key | Meaning |
|---|---|
case.tes_enabled |
Allow investment in thermal energy storage |
case.dh_connection |
Allow industry ↔ district heating heat exchange |
case.num_sectors |
Number of industry sectors (paper: 3; profiles exist for 1–3) |
case.investment_periods |
5 (paper) or 1. Only these two — the growth trajectories are five-element series |
case.co2_cap |
CO2 cap per period [t CO2], Figure 1 of the paper |
case.horizon_hours |
Truncate each period to this many hours (multiple of 60). 0 = full 8760. Smoke tests only — a shortened year distorts seasonal storage and annualised costs |
solver.name |
"HiGHS" or "Gurobi" |
output.csv_delim, output.csv_decimal |
Default ; and , (German Excel locale) |
run_case.jl Run one scenario
run_all_cases.jl Run the four paper scenarios
config.toml Paper configuration (3 sectors, 5 investment periods, full year)
config.quick.toml Smoke-test configuration (1 sector, 1 period, 30 days)
scripts/
explore_gui.jl Optional interactive visualisation of exported results
src/
IndustryDHFlex.jl Module entry point
inputs.jl Reading profiles and parameters from CSV
case.jl Resources, nodes, links, time structure
solve.jl Model creation and optimisation
results.jl Extracting and exporting solution values
gui.jl Optional EnergyModelsGUI visualisation
config.jl Reading config.toml
Data/
input_profiles.csv Hourly profiles, 26 series x 8760 h
parameters/ Investment, storage, OPEX, capacity and scalar parameters
design/top_level.yml Node positions for the optional GUI
results/ Model output (git-ignored)
All inputs are plain CSV. They were exported once from the original Excel workbooks and verified value-by-value against them.
Excel is deliberately not used at runtime: XLSX readers return cached formula values, so editing a parameter in a workbook without reopening it in Excel would silently feed stale numbers to the model. The original parameter sheet was also built on spilled dynamic-array formulas, which LibreOffice cannot recalculate at all. With CSV, what you read is what the model uses.
Data/input_profiles.csv — 8760 hourly values per series: electricity price, wind and
solar generation, gas and coal CHP production profiles, building heat demand, ambient
temperature, and per-sector demand ({DHHeat,LTHeat,HTHeat,Power}DemandSec1..3) and
excess heat (ExHeatSec1..3, DHExHeatSec1..3).
Data/parameters/ — investments.csv (capex, lifetime), storage_investments.csv,
opex.csv (fixed and variable), starting_capacities.csv, scalars.csv (discount rate,
gas price, emission price, transmission cost).
The profiles were derived from public data. The intermediate spreadsheets used to build them are not redistributed here; the underlying sources are:
| Input | Source |
|---|---|
| Hourly electricity price, wind/solar generation, gas & coal CHP profiles | Energy-Charts, Germany 2024 |
| District heating production shares and renewable share | BDEW Fernwärme statistics |
| Storage cost and performance parameters | Danish Energy Agency Technology Catalogue — Energy Storage |
| Industrial excess heat potentials | sEEnergies excess heat database |
| Grid mix emission factor (0.38 tCO2/MWh) | Umweltbundesamt 23/2024 |
| Gas (0.200) and hard coal (0.338) emission factors | Quaschning, specific CO2 emissions dataset |
Industrial heat demand profiles are described in Section 3 of the paper.
Five things in the case setup are easy to misread from the code alone.
Energy in MWh, power and capacity in MW, all costs in kEUR. So GasPrice = 0.0734
is 73.4 EUR/MWh, and a capex of 89.58 is 89.58 kEUR/MW. Note that kEUR/MW and EUR/kW
are numerically identical, which is why the paper's tables and the CSVs agree despite
using different unit labels.
The resource names in resources() are offset by one level from the temperature labels
used in the paper. The storage node names do follow the paper. Mapping:
| Code resource | Temperature | Paper level |
|---|---|---|
ExHeat |
40 °C | Excess heat |
DHHeat |
80 °C | Low temperature (LT) |
LTHeat |
150 °C | Medium temperature (MT) |
HTHeat |
300 °C | High temperature (HT) |
So LTHeatStorage stores DHHeat and is the paper's LT TES / tank-TES;
MTHeatStorage stores LTHeat and is the paper's MT TES (steam accumulator);
HTHeatStorage stores HTHeat and is the paper's HT TES (molten salt). The paper
gives the high-temperature level as ">200 °C"; the model uses 300 °C, which is what the
heat-pump COP arithmetic is based on.
Each industry sector draws electricity through a CapacityCostLink
(EnergyModelsFlex), not a plain Direct link. It divides the year into 60 sub-periods
of ~146 h, takes the maximum hourly import in each, and charges cap_cost
(1 kEUR/MW) on it — up to 60 kEUR/MW per year if the peak is flat across sub-periods.
This is a demand charge in the region of a German industrial Leistungspreis, and it is
one of the drivers of the peak-shaving behaviour reported in Section 3.1 of the paper.
District heating is connected to TransmissionConstraint by a plain link and pays no
capacity charge. The 93 EUR/MWh energy add-on (TransmissionCost) applies to both and
is levied at the TransmissionConstraint node.
Version-sensitive. The third numeric argument (
60) means number of sub-periods in EnergyModelsFlex 0.2.x. A later breaking release renamed the field and changed its meaning to duration, so the same60would silently become a 60-hour period. The[compat]bound inProject.tomlpins 0.2.x for this reason — do not relax it without reading the EnergyModelsFlex migration notes.
All RefSink demand nodes carry :surplus and :deficit penalties of 100 — that is
100 kEUR/MWh, three orders of magnitude above the most expensive supply option. The
slacks exist so that an infeasible configuration fails visibly rather than silently, and
are never used in the reported runs. Note that they are not among the variables exported
to combined_*.csv; set output.all_variables = true to inspect them.
| Quantity | Limit | Where |
|---|---|---|
Electricity through TransmissionConstraint |
500 MW, fixed, no investment option | single_nodes |
| Industry–DH heat exchange per sector | 100 MW total | investment["DHConnection"] |
Pit TES (PTESDH) level |
40 000 MWh total, 10 000 MWh added per period | investment_storage["SeasonalTESDH"] |
| Tank TES and industry LT TES level | 3 000 MWh added per period, no total cap | investment_storage["LTHeatStorage"] |
| Wind, solar | 300 MW added per period | investment["Wind"], ["Solar"] |
The CO2 cap is annual (emissions in one modelled year), not a five-year budget. The first period's cap of 1e9 t is deliberately non-binding and serves as an unconstrained reference.
Each run writes to results/<case name>/:
| File | Contents |
|---|---|
combined_op_df.csv |
Operational results, wide: one column per variable × node, one row per hour. First row is the column sum. |
combined_st_df.csv |
Strategic (per investment period) results, same layout |
metadata_op.csv, metadata_st.csv |
Case name, variable and node behind each column of the corresponding combined file |
all_variables.csv |
Complete dump of every model variable. Only when output.all_variables = true — see below |
The combined files use ; as separator and , as decimal mark (German Excel locale);
the metadata files use plain ,. Change this with output.csv_delim and
output.csv_decimal.
Because the combined files have no header row, metadata_*.csv is how you identify
columns: row i of the metadata names column i of the combined file.
Set output.per_variable_files = true to also get one CSV per variable under
results/<case name>/Cases/<case name>/.
The combined_*.csv files are selective, and their per-node-type dispatch does not
cover FlexibleOutput. In practice that means four node types per sector export no
columns at all:
GasBoilerSec<n>, CCSSec<n>, ElectricBoilerSec<n>, TempDropSec<n>
These nodes are fully present in the optimisation and constrain the solution normally —
they are simply absent from the exported tables. Their district-heating counterparts
(GasBoilerDH, ElectricBoilerDH) are RefNetworkNode and do export. Consequently
aggregates such as el_consumption and LT_heat_production computed from
combined_op_df.csv exclude the industrial boilers.
This is inherited behaviour, preserved deliberately so the exported files match those
used for the paper. If your analysis needs those nodes, set output.all_variables = true
to get a complete, unfiltered dump of every model variable:
variable,index,value
cap_use,n_ElectricBoilerSec1|sp1-t1,0.0
It is large — roughly 15 MB for the quick config, and much larger for the full horizon.
julia --project -e 'using Pkg; Pkg.add("EnergyModelsGUI")'
julia --project scripts/explore_gui.jlEnergyModelsGUI is not a dependency because it pulls in GLMakie/GeoMakie and needs a
display; excluding it keeps the default install roughly three times smaller. Node
positions live in design/top_level.yml.
- The published numbers have not been re-derived for this release. The code is verified to build, solve and export, and to be numerically identical to the version that produced the manuscript — but the headline results were generated with Gurobi, and that has not been re-run here. Treat the paper as the authority on the numbers.
investment_periodsaccepts only 1 or 5. The demand growth, renewable capacity and battery cost trajectories are hard-coded five-element series (Section 3 of the paper); other horizons would need those extended.- The
combined_*.csvexports omitFlexibleOutputnodes — seeall_variables.csv. The optimisation is unaffected; only the exported tables are. config.quick.tomltruncates the year to 30 days. It verifies the machinery, not the physics: seasonal storage and annualised costs are meaningless at that horizon.Project.tomlcompat bounds are load-bearing, not hygiene. At least one argument ofCapacityCostLinkchanged meaning in a breaking release of EnergyModelsFlex after the pinned 0.2.x — see Capacity charge on industrial electricity. Relaxing the bounds can change the model without changing this repository.- Perfect foresight, a single social-planner perspective, and linearised physics — the standard limitations of this model class, discussed in Section 2 of the paper.
If you use this code, please cite both the paper and the software. See
CITATION.cff, or use GitHub's "Cite this repository" button.
MIT — see LICENSE. Copyright (c) 2026 SINTEF Energy Research.