Energy System Optimization with Julia
Hamburg University of Applied Science - Summer 2025
Objective: Minimize total annual costs (investment + operational) while respecting budget constraints for multi-energy systems
Decision Variables:
Key Constraints:
Note
The Multi-Energy System Design problem extends the basic energy system design by incorporating multiple energy carriers (electricity and hydrogen) and energy conversion technologies (electrolyzers). The stochastic formulation considers multiple scenarios to account for uncertainty in renewable generation, demand patterns, and price variations across different energy carriers.
\(\text{Minimize} \quad \sum_{s \in \mathcal{S}} AC^{inv}_s + \sum_{w \in \mathcal{W}} AC^{inv}_w + \sum_{v \in \mathcal{V}} AC^{inv}_v + \sum_{x \in \mathcal{X}} AC^{inv}_x + \sum_{\omega \in \Omega} \pi_{\omega} (AC^{grid,imp}_{\omega} - AR^{grid,exp}_{\omega})\)
\(\sum_{s \in \mathcal{S}} (C^{E}_s e^{nom}_s + C^{P,ch}_s p^{ch,nom}_s + C^{P,dis}_s p^{dis,nom}_s) + \sum_{w \in \mathcal{W}} C^{W}_w p^{nom}_w + \sum_{v \in \mathcal{V}} C^{PV}_v p^{nom}_v + \sum_{x \in \mathcal{X}} C^{X}_x p^{nom}_x \leq B^{max}\)
\(\sum_{w \in \mathcal{W}} p_{w,t,\omega} + \sum_{v \in \mathcal{V}} p_{v,t,\omega} + (p^{in}_{t,\omega} - p^{out}_{t,\omega}) + \sum_{s \in \mathcal{S}_{\text{El}}} (p^{dis}_{s,t,\omega} - p^{ch}_{s,t,\omega}) = \sum_{x \in \mathcal{X}} p_{x,t,\omega} \quad \forall t \in \mathcal{T}, \omega \in \Omega\)
\(\sum_{x \in \mathcal{X}} h_{x,t,\omega} + \sum_{s \in \mathcal{S}_{\text{H2}}} (p^{dis}_{s,t,\omega} - p^{ch}_{s,t,\omega}) = d_{t,\omega} \quad \forall t \in \mathcal{T}, \omega \in \Omega\)
\(0 \leq p_{w,t,\omega} \leq f_{w,t,\omega} p^{nom}_w \quad \forall w \in \mathcal{W}, t \in \mathcal{T}, \omega \in \Omega\) \(0 \leq p_{v,t,\omega} \leq f_{v,t,\omega} p^{nom}_v \quad \forall v \in \mathcal{V}, t \in \mathcal{T}, \omega \in \Omega\) \(0 \leq p^{ch}_{s,t,\omega} \leq p^{ch,nom}_s \quad \forall s \in \mathcal{S}, t \in \mathcal{T}, \omega \in \Omega\) \(0 \leq p^{dis}_{s,t,\omega} \leq p^{dis,nom}_s \quad \forall s \in \mathcal{S}, t \in \mathcal{T}, \omega \in \Omega\) \(DoD_s e^{nom}_s \leq e_{s,t,\omega} \leq e^{nom}_s \quad \forall s \in \mathcal{S}, t \in \mathcal{T}, \omega \in \Omega\)
\(e_{s,t,\omega} = (1-sdr_s)e_{s,t-1,\omega} + \eta^{ch}_s p^{ch}_{s,t,\omega} - \frac{p^{dis}_{s,t,\omega}}{\eta^{dis}_s} \quad \forall s \in \mathcal{S}, t \in \mathcal{T}, \omega \in \Omega\)
\(0 \leq p_{x,t,\omega} \leq p^{nom}_x \quad \forall x \in \mathcal{X}, t \in \mathcal{T}, \omega \in \Omega\) \(h_{x,t,\omega} = \beta_x p_{x,t,\omega} \quad \forall x \in \mathcal{X}, t \in \mathcal{T}, \omega \in \Omega\)
Tip
The multi-energy system model can be adapted to specific use cases, i.e., if certain components already have fixed capacities, the corresponding variables can be converted to parameters. The stochastic formulation allows for more robust investment decisions by considering multiple scenarios across different energy carriers.
Note
These extensions enable comprehensive optimization of energy value chains, allowing for optimal sizing and operation of systems that convert between the different energy carriers while considering uncertainty and operational constraints.
Tip
You can ask questions anytime in class or via email!
Big M constraints are a fundamental modeling technique in Mixed-Integer Linear Programming (MILP) used to create conditional relationships between continuous and binary variables. The “M” represents a large constant that acts as an upper bound, allowing constraints to be “turned on” or “turned off” based on binary variable values.
Note
Big M constraints enable the modeling of complex logical relationships and conditional constraints that would otherwise be difficult to express in linear programming formulations.
The general form of a big M constraint is:
\(f(x) \leq M (1 - y)\)
where:
Logic: When \(y = 1\), the constraint becomes \(f(x) \leq 0\) (active). When \(y = 0\), the constraint becomes \(f(x) \leq M\) (inactive, since \(M\) is large).
Consider the constraint from the intermission lecture where generator 3 and generator 4 cannot start up within 3 time periods of each other:
\(\sum_{\tau \in [t, \min(t+3,|\mathcal{T}|)]} v_{3,\tau} \leq (1 - v_{4,t}) M \quad \forall t \in \mathcal{T}\)
\(\sum_{\tau \in [t-1, \max(t-3,1)]} v_{3,\tau} \leq (1 - v_{4,t}) M \quad \forall t \in \mathcal{T}\)
Tip
The choice of \(M\) is crucial: it must be large enough to make constraints inactive when needed, but not so large that it causes numerical issues in the solver.
For the startup timing constraint, a good choice for \(M\) would be:
Indicator constraints replace big M constraints with direct logical relationships. They are best suited for simple binary conditions where you have clear “if-then” relationships.
Big M approach:
\(p_g \leq p_g^{max} u_g\)
Indicator constraint approach:
@constraint(model, u_g => {p_g <= p_g_max}) # If generator is ON, then power ≤ max
@constraint(model, !u_g => {p_g == 0}) # If generator is OFF, then power = 0When to use: Simple binary logic (on/off, yes/no, active/inactive)
Key characteristic: Direct logical statements without mutual exclusivity requirements
Advantages:
Disadvantages:
Disjunctive programming handles complex multi-state conditions where exactly one state must be active. It uses special ordered sets (SOS1) to enforce mutual exclusivity.
Example: A generator can operate in exactly one of three modes: off, minimum load, or full load.
Big M approach:
\(p_g \leq (1 - u_g^{off}) + p_g^{min} u_g^{min} + p_g^{max} u_g^{max}\)
\(p_g \geq (1 - u_g^{off}) + p_g^{min} u_g^{min} + p_g^{max} u_g^{max}\)
\(u_g^{off} + u_g^{min} + u_g^{max} \leq 1\)
Disjunctive approach:
# SOS1 ensures exactly one binary variable is 1, others are 0
@constraint(model, [u_g_off, u_g_min, u_g_max] in MOI.SOS1(3))
@constraint(model, u_g_off => {p_g == 0}) # State 1: Off
@constraint(model, u_g_min => {p_g == p_g_min}) # State 2: Minimum load
@constraint(model, u_g_max => {p_g == p_g_max}) # State 3: Full loadWhen to use: Multiple mutually exclusive states or operational modes Key characteristic: Enforces that exactly one state is active at a time
Advantages:
Disadvantages:
| Aspect | Indicator Constraints | Disjunctive Programming |
|---|---|---|
| Complexity | Simple binary logic | Multiple states |
| Mutual Exclusivity | Not enforced | Automatically enforced |
| Use Case | “If A, then B” | “Must be in exactly one of states A, B, C” |
| Implementation | Direct logical statements | SOS1 constraints |
| Example | Generator on/off | Generator off/min/max modes |
Sometimes big M constraints can be reformulated using additional variables or constraints that are more efficient or numerically stable.
Example: Instead of using big M for minimum up time constraints, use a different approach.
Big M approach:
\(\sum_{t'=t}^{t+T^{min}-1} u_{g,t'} \geq T^{min} v_{g,t}\)
Reformulation approach:
# Use a continuous variable to track remaining minimum up time
@variable(model, 0 <= remaining_up_time[g, t] <= T_min)
@constraint(model, remaining_up_time[g, t] <= T_min * u_g[t])
@constraint(model, remaining_up_time[g, t] <= remaining_up_time[g, t-1] - 1 + T_min * v_g[t])
@constraint(model, remaining_up_time[g, t] >= 0)Advantages:
Disadvantages:
Note
Big M constraints are a powerful tool in MILP modeling, but they should be used carefully with appropriate values for \(M\) to ensure both correctness and computational efficiency.
A multi-energy system is a system that combines multiple energy carriers (e.g. electricity, hydrogen, heat, etc.) to meet the energy demand of a specific application.
Example case of a hydrogen value chain
Note
The figure shows a multi-energy system with:
In this lecture, we reformulate the Stochastic Multi-Energy System Design problem to be applicable to a wide range of energy system applications only by changing the sets and parameters. This allows us to:
Note
Remember: The design problem combines investment planning (sizing) with operational optimization, enabling comprehensive (multi-)energy system design.
… extending the stochastic design problem to a general multi-energy system framework
Note
The general framework uses technology-agnostic sets that can represent any energy conversion technology (electrolyzers, fuel cells, heat pumps, etc.) and any energy carrier (electricity, hydrogen, heat, natural gas, etc.).
Maximize the following components:
\(\text{Maximize}\)
Revenue from selling energy carriers:
\(\sum_{\omega \in \Omega} \pi_{\omega}^{prob} \sum_{e \in \mathcal{E}} \sum_{t \in \mathcal{T}} \pi_{e,t,\omega}^{sell} f_{e,t,\omega}^{sell}\)
Revenue from satisfied demand:
\(+\sum_{\omega \in \Omega} \pi_{\omega}^{prob} \sum_{g \in \mathcal{G}^{Dem}} \sum_{t \in \mathcal{T}} \pi_g^{dem} D_{e_g,t,\omega}\)
Costs of buying energy carriers:
\(-\sum_{\omega \in \Omega} \pi_{\omega}^{prob} \sum_{e \in \mathcal{E}} \sum_{t \in \mathcal{T}} \pi_{e,t,\omega}^{buy} f_{e,t,\omega}^{buy}\)
Taxes and levies costs:
\(-\sum_{\omega \in \Omega} \pi_{\omega}^{prob} \sum_{g \in \mathcal{G}^{Sup}} \sum_{e \in \mathcal{E}} \sum_{t \in \mathcal{T}} \tau_g f_{e,t,\omega}^{buy} \quad \text{where } e_g^{sup} = e\)
Variable operating costs:
\(-\sum_{\omega \in \Omega} \pi_{\omega}^{prob} \sum_{g \in \mathcal{G}} \sum_{t \in \mathcal{T}} C_g^{var} x_{g,t,\omega}^{total}\)
Startup costs:
\(-\sum_{\omega \in \Omega} \pi_{\omega}^{prob} \sum_{g \in \mathcal{G}^{UC}} \sum_{t \in \mathcal{T}} c_{g,t,\omega}^{start}\)
Investment costs:
\(-\sum_{g \in \mathcal{G}} AC_g^{inv}\)
Note
The objective maximizes profit: revenue from selling energy carriers and satisfied demand minus costs of buying energy carriers (including taxes and levies), variable operating costs, startup costs, and investment costs.
\(AC_g^{inv} = \frac{C_g^{inv}}{F^{PVAF}} P_g^{nom} \quad \forall g \in \mathcal{G} \setminus \mathcal{G}^S\)
\(AC_g^{inv} = \frac{C_g^{inv}}{F^{PVAF}} P_g^{nom} + \frac{C_g^{SOC}}{F^{PVAF}} SOC_g^{nom} \quad \forall g \in \mathcal{G}^S\)
\(\sum_{g \in \mathcal{G}} AC_g^{inv} \leq B^{max}\)
\(x_{g,e,t,\omega}^{in} = \sigma_{g,e}^{in} x_{g,t,\omega}^{total} \quad \forall g \in \mathcal{G}, e \in \mathcal{E}_g^{In}, t \in \mathcal{T}, \omega \in \Omega\)
\(x_{g,e,t,\omega}^{out} = \sigma_{g,e}^{out} \Theta_g x_{g,t,\omega}^{total} \quad \forall g \in \mathcal{G} \setminus \mathcal{G}^S, e \in \mathcal{E}_g^{Out}, t \in \mathcal{T}, \omega \in \Omega\)
\(x_{g,t,\omega}^{total} \leq P_g^{nom} \quad \forall g \in \mathcal{G} \setminus \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
\(x_{g,t,\omega}^{total} \leq M u_{g,t,\omega} \quad \forall g \in \mathcal{G} \setminus \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
\(0.1 P_g^{nom} \leq x_{g,t,\omega}^{total} + M (1 - u_{g,t,\omega}) \quad \forall g \in \mathcal{G} \setminus \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
where \(M\) is a large constant (big M).
Note
The minimum power constraint is implemented using big M constraints to maintain linearity. When \(u_{g,t,\omega} = 1\) (unit is on), the third constraint becomes \(0.1 P_g^{nom} \leq x_{g,t,\omega}^{total}\), enforcing the minimum power requirement. When \(u_{g,t,\omega} = 0\) (unit is off), the second constraint forces \(x_{g,t,\omega}^{total} = 0\).
\(0 \leq x_{g,t,\omega}^{total} \leq f_{g,t,\omega} P_g^{nom} \quad \forall g \in \mathcal{G}^{Var}, t \in \mathcal{T}, \omega \in \Omega\)
where \(f_{g,t,\omega}\) is the capacity factor for variable technology \(g\).
\(soc_{g,t,\omega} = (1 - \delta_g) soc_{g,t-1,\omega} + \eta_g^{ch} x_{g,t,\omega}^{total} - \frac{\sum_{e \in \mathcal{E}_g^{Out}} \sigma_{g,e}^{out} x_{g,e,t,\omega}^{out}}{\eta_g^{dis}} \quad \forall g \in \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
\(0 \leq soc_{g,t,\omega} \leq SOC_g^{max} \quad \forall g \in \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
\(soc_{g,t=0,\omega} = soc_{g,t=|\mathcal{T}|,\omega} = SOC_g^{init} \quad \forall g \in \mathcal{G}^S, \omega \in \Omega\)
\(x_{g,t,\omega}^{total} \leq P_g^{nom} \quad \forall g \in \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
\(x_{g,t,\omega}^{total} \leq M o_{g,t,\omega} \quad \forall g \in \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
\(\sum_{e \in \mathcal{E}_g^{Out}} \sigma_{g,e}^{out} x_{g,e,t,\omega}^{out} \leq P_g^{nom} \quad \forall g \in \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
\(\sum_{e \in \mathcal{E}_g^{Out}} \sigma_{g,e}^{out} x_{g,e,t,\omega}^{out} \leq M (1-o_{g,t,\omega}) \quad \forall g \in \mathcal{G}^S, t \in \mathcal{T}, \omega \in \Omega\)
where \(M\) is a large constant (big M).
Note
The storage charging/discharging mutual exclusion is implemented using big M constraints. When \(o_{g,t,\omega} = 1\) (charging mode), the second constraint allows charging while the fourth constraint forces discharging to zero. When \(o_{g,t,\omega} = 0\) (discharging mode), the second constraint forces charging to zero while the fourth constraint allows discharging.
\(x_{g,e_g^{sup},t,\omega}^{out} = Q_g^{nom} \quad \forall g \in \mathcal{G}^{Sup}: \text{supply type} = \text{"fix"}, t \in \mathcal{T}, \omega \in \Omega\)
\(x_{g,e_g^{sup},t,\omega}^{out} \leq Q_g^{nom} \quad \forall g \in \mathcal{G}^{Sup}: \text{supply type} = \text{"max"}, t \in \mathcal{T}, \omega \in \Omega\)
\(\sum_{g \in \mathcal{G}} x_{g,e,t,\omega}^{out} + f_{e,t,\omega}^{buy} = \sum_{g \in \mathcal{G}} x_{g,e,t,\omega}^{in} + f_{e,t,\omega}^{sell} + D_{e,t,\omega} \quad \forall e \in \mathcal{E}^D, t \in \mathcal{T}, \omega \in \Omega\)
\(\sum_{g \in \mathcal{G}} x_{g,e,t,\omega}^{out} + f_{e,t,\omega}^{buy} = \sum_{g \in \mathcal{G}} x_{g,e,t,\omega}^{in} + f_{e,t,\omega}^{sell} \quad \forall e \in \mathcal{E} \setminus \mathcal{E}^D, t \in \mathcal{T}, \omega \in \Omega\)
\(v_{g,t,\omega} \geq u_{g,t,\omega} - u_{g,t-1,\omega} \quad \forall g \in \mathcal{G}^{UC}, t \in \mathcal{T}, \omega \in \Omega\)
\(c_{g,t,\omega}^{start} \geq C_g^{start} v_{g,t,\omega} \quad \forall g \in \mathcal{G}^{UC}, t \in \mathcal{T}, \omega \in \Omega\)
Note
This general framework can be specialized to specific applications by:
Questions?
For more interesting literature to learn more about Julia, take a look at the literature list of this course.
Lecture XI - Framework for Multi-Energy System Optimization | Dr. Tobias Cors | Home