Energy System Optimization Models

Updated 9 February 2026

Energy System Optimization Models (ESOMs) are advanced mathematical frameworks that optimize capacity expansion, dispatch, sector coupling, and infrastructure investments in multi-energy systems.
They employ linear and mixed-integer programming formulations to incorporate detailed temporal, spatial, and discrete constraints that enable rigorous decarbonization planning.
Modern ESOM frameworks integrate sector coupling, uncertainty modeling, and decomposition techniques to efficiently solve large-scale, multi-period optimization problems.

Energy system optimization models (ESOMs) are mathematical models—typically large-scale linear programs (LPs) or mixed-integer linear programs (MILPs)—that provide rigorous frameworks for optimizing capacity expansion, operational dispatch, sectoral integration, and infrastructure development in complex multi-energy systems. ESOMs are foundational in long-term decarbonization planning, political decision support, and techno-economic evaluation of future energy scenarios. They enable system-wide co-optimization of investments and operations, explicitly accounting for detailed techno-economic, temporal, and spatial constraints, and have evolved to handle increasing model sizes, computational complexity, sector coupling, and uncertainty (Miehling et al., 2023, Göke, 2020).

1. Mathematical Structure and Core Formulation

The canonical ESOM is a multi-period optimization that minimizes total system cost, subject to energy balances, capacity constraints, storage dynamics, unit-commitment, policy constraints (e.g., emissions caps), and potentially network flow limits. The decision vector includes variables for technology investments, dispatch, storage operations, start-up/shutdown events, and network expansions. The standard LP/MILP representation is

$\min_{x} \; c^T x \qquad \text{s.t.} \;\; Ax \leq b,\quad x \geq 0,\quad x_j \in \mathbb{Z} \;\; \text{for discrete variables.}$

Core constraints include:

Carrier balances (electricity, heat, fuels) at each node and timestep.
Technology- and time-dependent upper and lower dispatch (e.g., $0 \leq g_{i,n,t} \leq \rho_{i,n,t} \kappa_{i,p(t)}$ ).
Storage state-of-charge recursions.
Unit-commitment binaries with minimum up/down times, ramping, and start-up/shutdown coupling.
Emissions caps, renewable portfolio standards, and multi-period capacity evolution (vintaging/depreciation).
Investment and flow bounds, often referencing spatial and temporal aggregations (Miehling et al., 2023, Göke, 2020, Schönfeldt et al., 2020, Göke, 2020, Schönfeldt et al., 2020).

Colloquially, the structural block-matrix layout (arrowhead, staircase, or block-diagonal) is exploited for computational performance (Bröchin et al., 2022, Hadidi et al., 29 Jul 2025).

2. Temporal and Spatial Aggregation

Temporal disaggregation to hourly resolution is computationally prohibitive in continental or multi-sectoral studies; current practice is to aggregate time series into representative periods (RPs) using clustering methods (k-means, hierarchical, or model-adaptive) (Zhang et al., 2022). Traditional “block-diagonal” aggregation fails to respect chronology at RP boundaries, breaking continuity for storage and commitment variables. The Markov-matrix linkage method (Auer et al., 21 Oct 2025) addresses this by introducing a representative-period transition probability matrix $P$ :

State space $S = \{1,\ldots, RP\}$ ; empirically count $n_{ij}$ for consecutive period transitions.
Transition matrix $P_{ij} = n_{ij} / \sum_{k} n_{ik}$ .
Chronological boundary variables (e.g., storage $l_{rp,1}$ , unit on/off status) are linked via expectations: $\mathbb{E}[x_{rp,1}] = \sum_{i=1}^{RP} P_{i,rp} x_{i,K}$ .
For discrete variables at RP edges, integer constraints are relaxed (fractional expectation), preserving feasibility and keeping errors to <1% in commitment metrics—an order of magnitude improvement over “cyclic” linking, with minimal runtime overhead (Auer et al., 21 Oct 2025). This approach enables highly aggregated, yet chronologically accurate, tractable ESOMs.

Spatial aggregation employs clustering on demand, resource, or network features (e.g., k-means, spectral clustering, impedance-based), balancing model size and fidelity. Adaptive clustering using ESOM-native decision variable sensitivities further reduces error in approximated solutions (Zhang et al., 2022).

3. Model Classes and Framework Implementations

Modern ESOM frameworks (e.g., AnyMOD.jl, PyPSA-Eur, ETHOS.FINE, FOCUS, OpTUMus) implement graph- or data-frame–backed architectures for modularity, multi-period expansion, sectoral coupling, and multi-carrier representation:

AnyMOD.jl uses a directed graph abstraction for hierarchical index sets (multi-level time/carrier/region), native JuMP embedding, and per-carrier or per-node temporal/spatial granularity (Göke, 2020).
PyPSA(-Eur, -Earth) establishes explicit network representation, time-and-spatial clustering, and sector-coupled optimization in a fully open-source stack, using linear programming by default (Parzen et al., 2022).
ETHOS.FINE supports integrated multi-carrier, multi-region, multi-vector pathways with full flexibility in space, time, and investment epoch scaling (Groß et al., 2023).
FOCUS targets multi-layered (prosumer-district-city) urban energy modeling with MILP, multi-objective front computation, and modular scenario configuration (Gong et al., 2023).
OpTUMus offers a parameterized constraint-builder abstraction, separating data layers from symbolic LP/MILP structure for technological and sectoral extensibility (Miehling et al., 2023).

Common features include CSV- or pandas-based data ingestion, open-source solver integration (Gurobi, CPLEX, HiGHS, GLPK, CBC), and built-in visualization and post-processing tools.

4. Computational Complexity and Solution Strategies

The scale and structure of ESOMs—driven by spatiotemporal granularity, discrete variables, and scenario expansion—produce models with tens to hundreds of millions of variables/constraints (Kempke et al., 26 Feb 2025, Sasanpour et al., 14 Aug 2025, Hadidi et al., 29 Jul 2025). Three main approaches optimize tractability:

Numerical conditioning and scaling: Automatic scaling per physical unit (e.g., power-of-two shifts by unit type) dramatically reduces solver time by 2–3× by controlling the coefficient range in the constraint matrix (from $\sim 10^{16}$ to $\sim 10^{9}$ ), directly mitigating roundoff and simplex/barrier stalling (Bröchin et al., 2022, Göke, 2020).
Decomposition and parallelization:
- Benders decomposition, Dantzig–Wolfe, and Lagrangian relaxation exploit block structure, enabling decomposition by scenario, time, or region (Hadidi et al., 29 Jul 2025, Sasanpour et al., 14 Aug 2025).
- Temporally Split Benders Decomposition (TSBD) partitions subproblems along $T$ , complementing the classical scenario decomposition: each temporal block introduces storage linking variables, supporting massive parallelization (e.g., 60–80% runtime reduction on 87M variable German case studies) (Sasanpour et al., 14 Aug 2025).
Heuristics and LP-MIP relaxations: For extremely large-scale MILPs (e.g., 83M variables/900k binaries), heuristic techniques—relaxation-enforced neighbor search, ML-guided rounding, fix-and-propagate—can find feasible solutions 20–100× faster than state-of-the-art solvers at <1–35% primal-dual gap, often sufficient for scenario analysis (Kempke et al., 26 Feb 2025).

Solver choice, barrier-only vs. simplex + crossover modes, tolerance setting, threading, and model structure (sparsity, degeneracy, arrowhead block forms) all critically affect tractability (Bröchin et al., 2022).

5. Sector Coupling, Advanced Constraints, and Material-Energy Feedback

Leading ESOMs enable rigorous multi-sector modeling, including electricity, heat at multiple temperature levels, hydrogen, gas, carbon, and material flows:

Sectoral Coupling: Sector-integrated ESOMs use multi-carrier graphs (nodes for energy carriers, edges for conversion technologies) to enforce cross-carrier balances (e.g., electrolyzer/heat pump/cogeneration) (Göke, 2020). Temporal/hierarchical resolution can be mixed per sector (e.g., hourly for electricity, daily for heat).
Advanced technological constraints: Discrete-temperature linearization for heat demand, temperature-dependent source efficiency (heat pumps, STC), stratified thermal storage, and non-ideal dynamics are handled via discretization and optional binary assignment for tractable LPs (Schönfeldt et al., 2020).
Gas/power system coupling: Tight piecewise-linear MIP formulations for gas flow and linepack preserve solution quality while reducing binary variable count by >50%, enabling efficient integration of nonlinear gas dynamics in MILP ESOMs (average 2.8× speedup, <0.5% error in dispatch) (Klatzer et al., 2024).
Material demand feedback: Bulk-material needs for infrastructure (steel, cement, chemicals) are endogenous, creating an energy–material loop in cost-optimal planning. When material demand/production is integrated as a co-optimized block (e.g., PyPSA-Eur), investment shifts toward less material-intensive technologies, modifies hydrogen sourcing, and alters the mix of industrial routes—a phenomenon missed by models with exogenous material flows (Burghardt et al., 27 Aug 2025). Quantitatively, steel demand rises by +9% and utility PV falls up to 9% in fully endogenous scenarios.

6. Uncertainty Modeling, Learning, and Alternative Generation

Contemporary ESOMs account for uncertainty and innovation dynamics using multiple methodological advances:

Stochastic & robust optimization: Two-stage stochastic programming embeds scenario-based uncertainty; Benders decomposition is applied for scenario and time-splitting (Sasanpour et al., 14 Aug 2025). Robust and distributionally robust frameworks (with ambiguity sets or moment constraints) offer control over conservatism but increase computational complexity (Mabrouk et al., 5 Dec 2025).
Endogenous learning: Technological “learning-by-doing” is captured via experience curves (Wright-type, one- or multi-factor with per-component granularity). Nonlinearities are often handled with piecewise-linearization. Credible parameterization (learning rates, cost floors, build-rate caps) and tractable reformulations are vital for system-level realism (Ouassou et al., 2021).
Exploring near-optimal alternatives: Modeling to Generate Alternatives (MGA) is incorporated in mainstream ESOMs (e.g., Temoa) to systematically probe the near-optimal solution space under deep uncertainty. MGA formulations add a cost-slack constraint and iteratively seek maximally different (orthogonal) solutions via sequentially updated penalty objectives, revealing knife-edge trade-offs and robustness boundaries that single-solution models obscure (DeCarolis et al., 2019).

7. Model Application, Validation, and Limitations

ESOMs facilitate scenario-based decarbonization studies at national, continental, or global scales:

Systems such as PyPSA-Earth have established high-resolution, globally open ESOMs for multi-region or national planning, validated for technical demand, resource, and network data against official datasets (IRENA, TCN, GWEC) (Parzen et al., 2022).
Comparison of ESOM-approximate seasonal planning with full unit-commitment shows that, with appropriate constraint modifications, cost and dispatch errors remain <2% while yielding runtime reductions of ∼100× (Queiroz et al., 2019).
Limitations include data intensity (material and process parameters), assumption of steady-state or perfect foresight, challenges modeling nonconvexities, and the need for credible uncertainty quantification. Ongoing research addresses multi-period pathways, explicit recycling/retirement, stochastic dynamics, and scalable robust optimization (Burghardt et al., 27 Aug 2025, Mabrouk et al., 5 Dec 2025).

Key recent references:

Markov-matrix period linkage for RP-based MILPs (Auer et al., 21 Oct 2025)
Model-driven adaptive time clustering (Zhang et al., 2022)
Endogenous learning in capacity expansion (Ouassou et al., 2021)
Automated scaling and solver recommendations (Bröchin et al., 2022)
Gas–power MILP integration (Klatzer et al., 2024)
Temporally split Benders for stochastic ESOMs (Sasanpour et al., 14 Aug 2025)
Material–energy feedback integration (Burghardt et al., 27 Aug 2025)
Parallel decomposition survey (Hadidi et al., 29 Jul 2025)
Framework architectures (AnyMOD.jl (Göke, 2020), PyPSA-Earth (Parzen et al., 2022), ETHOS.FINE (Groß et al., 2023), FOCUS (Gong et al., 2023), OpTUMus (Miehling et al., 2023))