Deep Reinforcement Learning Approaches the MILP Optimum of a Multi-Energy Optimization in Energy Communities

1. Introduction

As modern energy systems become increasingly dominated by renewable, variable power generation, local energy communities (ECs) emerge as a promising way to stabilize the grid and utilize electricity where it is generated [1]. These advantages foster the transition from individual buildings toward Net Zero Energy Buildings [2] or even a Net Zero Energy Neighborhood [3,4]. For the optimization and ideal operation of such systems, Mixed-Integer Linear Programming (MILP) has been regarded as the gold standard for many years. MILP is a powerful method for finding global optima, making it widely used in model-based optimization across economics, logistics, and energy-related fields [5,6,7]. However, as the complexity of the investigated systems increases, the required linearizations become more intricate, and computation times grow exponentially due to the introduction of additional binary variables. This prompts the need for alternative techniques [8]. Recent advances in machine learning and computational power have renewed interest in alternative methods for optimization and control. Reinforcement learning (RL), being model-free and data-driven, significantly reduces computation time once an agent has been trained, while still being able to approach the optimal solutions found using MILP [9,10]. However, since reinforcement learning (RL) is data-driven, agents require a sufficient amount of high-quality data for training. Consequently, training time becomes a critical factor in real-world applications, especially as retraining may be necessary during deployment. And while model-based approaches do not depend on historical data, RL agents tend to outperform traditional methods in the optimization process due to their faster execution and greater computational efficiency once trained. In ECs, where distributed flexibilities can be aggregated and coordinated across multiple entities, RL-based optimization approaches present a promising solution. Due to their decentralized and dynamic nature, such systems benefit from RL agents’ ability to support real-time decision-making under uncertainty. Ideally, this leads to reduced costs for all participants [11]. Optimizing single households often requires transparency and a willingness to compromise on comfort, as optimal control may lead to small temperature variations or limited availability during unforeseen events. Neighborhood-level solutions can help mitigate these issues, as aggregation typically leads to a smoother load curve in larger systems [12,13,14]. In the following subsection, relevant publications on applications of MILP and RL optimization in different energy systems are critically discussed.

1.1. Related Works

In several studies on energy systems, MILP optimization has been investigated and compared to a rule-based conventional operation. Baumann et al. [15] used a co-simulation approach of IDA ICE and Gurobi to optimize the energy system’s control. The building was then optimized on a daily basis using MILP and simulated using a physical resistor-capacitor (RC) model. The investigation demonstrated significant benefits including reduced electricity costs, improved self-consumption and self-sufficiency. In the work of Aguilera et al. [16], MILP was used to control large-scale heat pumps in a simulation utilizing thermal demand forecasts as well as heat pump performance maps. This optimization achieved cost savings in the operation using the digital twin model. Kepplinger et al. [17] also applied the MILP approach in a simulation, incorporating forecasting and state estimation methods to optimally control a heating element in a domestic hot water heater. In their later work, they validated their results using an experimental real-world setup [18]. Cosic et al. [19] extended the application of MILP optimization from operational scheduling to investment planning. In a comprehensive framework, they optimized the sizing and placement of PV and storage systems in a real Austrian municipality, evaluating multiple tariff and market scenarios to demonstrate the model’s adaptability and precision. A 15% reduction in total energy costs and a 34% reduction in ${\mathrm{CO}}_2$ emissions was achieved.

RL has also been utilized as an optimizing method for energy systems, although primarily using incentive functions not related to economics. Bachseitz et al. [20] investigated various RL algorithms in comparison to a rule-based control strategy for managing the heat pump of a multi-family building. Their findings indicate that while the RL agent was able to maintain the required storage temperatures, it fell short of matching the rule-based approach in maximizing PV self-consumption. Lissa et al. [21] implemented an RL agent guided by a reward structure based on comfort levels and energy consumption to manage the energy system of a simulated single-family home. They achieved energy savings of 8 % as well as an increased use of renewable energy compared to a rule-based approach. Rohrer et al. [22] applied RL in a lab experiment to evaluate its feasibility for demand response. Using six months of real-world data for training, their approach achieved considerable energy savings, once more emphasizing the practical potential of RL in real demand response scenarios. Franzoso et al. [23] highlighted the use of reinforcement learning to optimize multi-energy systems integrating renewable technologies. Their study demonstrated improved energy management by reducing emissions as well as operational costs, showing the versatility of RL in energy applications.

Langer and Völling investigated a system comprising a single-family home with a heat pump, PV system, and battery electric storage system, optimized using an RL approach [24]. Their RL solution achieved a performance close to the results of a MILP model developed in their earlier work [25]. The evaluation focused on user comfort, grid feed-in, and overall energy usage, but did not account for economic incentives. While this demonstrates the feasibility of RL for residential energy management, the absence of cost considerations limits its applicability in scenarios where financial optimization is critical.

While MILP and RL have been widely studied for building-level energy optimization, few works compare both methods in a comprehensive manner and using real-world data from an existing neighborhood energy system. Furthermore, most RL studies focus on comfort or energy savings, often ignoring economic incentives driven by real-time electricity prices. Finally, the leveraging of price signals and the flexibility of storage, while accounting for PV production, remains insufficiently explored at the community scale. These gaps are addressed in the present study by directly comparing MILP and RL for economic optimization of an EC using historical real-time pricing and measured on-site data. The detailed contributions of this work are described in the following subsection.

1.2. Contribution

This study comprehensively compares RL and MILP side by side under closely matched conditions using both measured and synthetic data for an existing EC, including demand, real-time electricity prices and PV generation profiles. The optimization focuses on exploiting the flexibility of an existing thermal energy storage (TES) by shifting the heat pump operation toward periods of low electricity prices while avoiding high-price intervals, thereby reducing overall energy costs. To enable a thorough comparison of the optimization methods, a reference scenario without flexibility was used, where heat demand is met directly. The individual contributions of this study are:

• Side-by-side comparison of MILP and RL, for the economic optimization of an existing EC using real-world input data and real-time electricity prices.
• Development of an RL-based control strategy for economic demand response, leveraging price signals to shift grid usage toward low-price periods by optimizing the operation of thermal storage and PV flexibility, aiming to approach the MILP-derived optimum.
• Comprehensive benchmark of RL and MILP optimization in a direct comparison and compared to a no-flexibility reference scenario. Within the benchmark cost savings, operational strategies, and robustness under realistic conditions, including PV variability and demand patterns, are assessed.

Therefore, this study provides practical insights into the performance and trade-offs of RL and MILP in energy system optimization, paving the way towards future control strategies for cost-efficient operation of ECs.

2. Methods

In this section, the system is presented and the used methods are highlighted, from general equations over the RL implementation and the MILP.

2.1. System

A real-world energy system was investigated in this study. The small EC, as illustrated in Figure 1, comprises three single-family homes that share both electrical and thermal energy infrastructure. The electrical side includes a grid connection point and rooftop photovoltaic panels. The thermal energy system is centered around a sensible TES, providing flexibility to the two heat-supply components: a 12 kW geothermal heat pump (HP) serving as the primary heat source, and a 6 kW auxiliary resistance heating element (HE) to cover peak loads. The storage tank has a capacity of 1120 liters of water and operates between 35°C and 55°C ($T_{\mathrm{stor}}$), allowing for approximately 26 kWh of storage capacity ($E_{\mathrm{stor}}$).

The direction of the energy flows is indicated by arrows in Figure 1. All flows labeled with the letter P represent electrical power exchanged between system components. These include solar power generation ($P_{\mathrm{PV}}$), the general electrical load of the homes ($P_{\mathrm{load}}$), and the electrical power consumption of the heating components, the heat pump ($P_{\mathrm{HP}}$) and the heating element ($P_{\mathrm{HE}}$). All of these flows are unidirectional, with the exception of the grid connection ($P_{\mathrm{grid}}$), which can either source electricity from the grid ($P_{\mathrm{grid},\mathrm{pos}}$) or feed in surplus electricity generated by the PV system ($P_{\mathrm{grid},\mathrm{neg}}$).

Thermal energy flows, denoted by $\dot{Q}$, indicate the direction of heat transfer between components. The heat generated by the heat pump (${\dot{Q}}_{\mathrm{HP}}$) and the heating element (${\dot{Q}}_{\mathrm{HE}}$) is stored in the TES. While the majority of the stored energy is ultimately used to meet the EC’s heating demand ($\dot{Q}\mathrm{out}$), heat losses to the environment ($\dot{Q}\mathrm{loss}$) must also be taken into account. It is noteworthy that in the investigated real energy system, a bypass option to directly cover the heating demand ${\dot{Q}}_{\mathrm{out}}$ using the heat pump and/or heating element is not available. While the heat load ${\dot{Q}}_{\mathrm{out}}$ from the TES supplies the required energy for space heating, domestic hot water consumption is accounted for within the synthetically generated electric load profile $P_{\mathrm{load}}$.

2.2. Data

The data used for the simulation and optimization of the EC is categorized as real-world (measured) on-site data, historical data, or synthetic data, as summarized in Table 1.

Measured data were logged for a period of one year from October 19, 2022. The sensors recorded data at a one-minute resolution, which were averaged and stored every 15 minutes, resulting in 96 samples per day and a total of $35,040$ samples for the entire year. The synthetic data was generated with the same temporal resolution to ensure consistency across all datasets. Historical data was resampled to this resolution when necessary.

The measured temperature data exhibited considerable noise due to the interference from the heat pump operation especially when used for cooling the buildings in the summer period, resulting in large fluctuations of the measured ground probe temperature profile. To ensure smooth simulation and optimization, as well as realistic COP and power estimates for both optimization methods, the data was preprocessed prior to any algorithmic use. Namely, a rolling average filter for 96 values was applied to smoothen the signal, resulting in a synthetic signal derived from the measured signal.

For the RL optimization, the seasonal classification was considered by assigning binary values to each time step, where a value of 1 indicates the active season. Seasons were defined as shown in Table 2.

For validation and to ensure robustness of the algorithms, a train-test split was applied to the data, as seen in Figure 2. This split was designed to include all seasons and their combinations within both sets. Specifically, the data for each month were divided such that in total 260 of the days were used for training and the remaining 105 for testing. The exception was the initial month, October 2022, which was exclusively allocated to training, and the final month, October 2023, exclusively to testing. This approach ensures comprehensive seasonal coverage during both training and testing phases.

2.3. Physical Model

The EC was modeled through a series of equations to calculate the used energy and energy flows at each timestep. For this model, the following assumptions were made:

• The specific heat capacity of water is equal to $c_{\mathrm{water}}=4.18\mathrm{kJ}/(\mathrm{kg}·\mathrm{K})$, irrespective of the TES temperature.
• Spatial temperature variations in the TES are neglected, making it a single node model.
• The mass balance is always fulfilled for the TES. Due to the narrow operating temperature range of only 20 °C, the volume of the storage is assumed to remain constant, as thermal expansion effects are negligible.
• The electric power of the heating element is equivalent to its heating power ($\eta =1$).
• The COP of the heat pump is equal to $0.5·{\mathrm{COP}}_{\mathrm{Carnot}}$. This is in accordance with Walden and Pedulla [32].

The heat capacity $C_{\mathrm{stor}}$ of the sensible water storage was calculated using the volume of the tank $V_{\mathrm{stor}}$, the specific density of water ${\rho }_{\mathrm{water}}$, and the specific heat capacity of water $c_{\mathrm{water}}$:

$$C_{\mathrm{stor}}=V_{\mathrm{stor}}·{\rho }_{\mathrm{water}}·c_{\mathrm{water}}$$

(1)

Using the specific heat capacity and the temporal change of temperature $\mathrm{d}{T}_{\mathrm{stor}}\left(t\right)/\mathrm{d}t$ of the storage tank at each timestep, the change in stored energy can be calculated. This was done by accounting for the heat flows to and from the tank, namely the heat supplied by the heating element ${\dot{Q}}_{\mathrm{HE}}$, the heat pump ${\dot{Q}}_{\mathrm{HP}}$, the heating demand ${\dot{Q}}_{\mathrm{out}}$, and the heat losses ${\dot{Q}}_{\mathrm{loss}}$:

$$C_{\mathrm{stor}}·{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{stor}}\left(t\right)}{\mathrm{d}t}}={\dot{Q}}_{\mathrm{HE}}\left(t\right)+{\dot{Q}}_{\mathrm{HP}}\left(t\right)-{\dot{Q}}_{\mathrm{out}}\left(t\right)-{\dot{Q}}_{\mathrm{loss}}\left(t\right)$$

(2)

To calculate the heating power ${\dot{Q}}_{\mathrm{HE}}$ of the heating element, the electrical power input $P_{\mathrm{HE}}$ and the efficiency $\eta$ of the system are required:

$${\dot{Q}}_{\mathrm{HE}}\left(t\right)=\eta ·P_{\mathrm{HE}}\left(t\right)$$

(3)

The heating power of the heat pump, denoted as ${\dot{Q}}_{\mathrm{HP}}$, was computed based on its electrical input power $P_{\mathrm{HP}}$ and an estimated coefficient of performance (COP) at each timestep. To this end, the ideal Carnot COP was used to model the temperature dependency of the heat pump COP, which depends on the heat sink temperature $T_{\mathrm{ub}}$ and the ground probe temperature $T_{\mathrm{GP}}\left(t\right)$.

$${\mathrm{COP}}_{\mathrm{carnot}}\left(t\right)={\displaystyle \frac{T_{\mathrm{ub}}\left(t\right)}{T_{\mathrm{ub}}\left(t\right)-T_{\mathrm{GP}}\left(t\right)}}$$

(4)

$${\dot{Q}}_{\mathrm{HP}}\left(t\right)=0.5·{\mathrm{COP}}_{\mathrm{carnot}}\left(t\right)·P_{\mathrm{HP}}\left(t\right)$$

(5)

For simplification, the heating demand ${\dot{Q}}_{\mathrm{out}}\left(t\right)$ was defined as the total heat consumption of the three separate houses in the energy system, ${\dot{Q}}_{\mathrm{H}1}$, ${\dot{Q}}_{\mathrm{H}2}$, and ${\dot{Q}}_{\mathrm{H}3}$. These individual demands were aggregated into a single total demand:

$${\dot{Q}}_{\mathrm{out}}\left(t\right)={\dot{Q}}_{\mathrm{H}1}\left(t\right)+{\dot{Q}}_{\mathrm{H}2}\left(t\right)+{\dot{Q}}_{\mathrm{H}3}\left(t\right)$$

(6)

The heat losses of the storage tank ${\dot{Q}}_{\mathrm{loss}}\left(t\right)$ were calculated using the heat transfer coefficient h, the total surface area of the tank A, and the temperature difference between the storage temperature $T_{\mathrm{stor}}\left(t\right)$ and the ambient temperature $T_{\infty }$:

$${\dot{Q}}_{\mathrm{loss}}\left(t\right)=h·A·\left(T_{\mathrm{stor}}\left(t\right)-T_{\infty }\right)$$

(7)

Using the previously defined heat capacity $C_{\mathrm{stor}}$ (see Equation (1)), the temporal rate of change of the storage tank temperature can also be expressed by:

$${\displaystyle \frac{{dT}_{\mathrm{stor}}\left(t\right)}{dt}}={\displaystyle \frac{{\dot{Q}}_{\mathrm{HE}}\left(t\right)+{\dot{Q}}_{\mathrm{HP}}\left(t\right)-{\dot{Q}}_{\mathrm{out}}\left(t\right)-{\dot{Q}}_{\mathrm{loss}}\left(t\right)}{C_{\mathrm{stor}}}}$$

(8)

Alternatively, assuming constant heat fluxes and parameters over time, the storage temperature evolution can be analytically described by the solution of the first-order linear differential equation:

$$\begin{array}{cc}\hfill T_{\mathrm{stor}}\left(t\right)& ={\displaystyle \frac{b}{a}}+\left(T_0-{\displaystyle \frac{b}{a}}\right)·e^{(-a·t)}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill a& ={\displaystyle \frac{hA}{C_{\mathrm{stor}}}}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill b& ={\displaystyle \frac{{\dot{Q}}_{\mathrm{HE}}+{\dot{Q}}_{\mathrm{HP}}-{\dot{Q}}_{\mathrm{out}}+hA·T_{\infty }}{C_{\mathrm{stor}}}}\hfill \end{array}$$

(11)

The heating power ${\dot{Q}}_{\mathrm{heat}}$ was controlled using a simple, discrete-time proportional P-controller with saturation, as shown in Algo. A1. At each time step i, with a duration of $\Delta t_{\mathrm{PI}}={\displaystyle \frac{1}{60}}\mathrm{h}$, the control error $e_i$ was computed as the difference between the desired storage temperature $T_{\mathrm{set},i}$ and the actual temperature $T_{\mathrm{actual},i}$. The heating power was then calculated by multiplying this error with the proportional gain $B_0$. To ensure that the power remains within operational limits, ${\dot{Q}}_{\mathrm{heat},i}$ was constrained to lie within the range defined by ${\dot{Q}}_{\mathrm{HP},\min }$ and ${\dot{Q}}_{\mathrm{HP},\max }$. The heating element was set to activate if the temperature of the TES drops to 2 degrees below $T_{\mathrm{lb}}$, acting strictly as an auxiliary heating unit. Both the RL and the MILP optimization are able to freely adjust the set temperature in order to utilize the flexibility of the thermal storage system. This flexibility is determined by the thermal capacity $C_{\mathrm{stor}}$ and the temperature bounds $T_{\mathrm{lb}}$ and $T_{\mathrm{ub}}$, set to 35°C and 55°C, respectively.

A state of charge (SOC) was defined for the TES, and it is computed based on the storage temperature at time t according to:

$$\mathrm{SOC}\left(t\right)={\displaystyle \frac{T_{\mathrm{ub}}-T_{\mathrm{stor}}\left(t\right)}{T_{\mathrm{ub}}-T_{\mathrm{lb}}}}$$

(12)

where $T_{\mathrm{stor}}\left(t\right)$ is the current storage temperature, and $T_{\mathrm{ub}}$ and $T_{\mathrm{lb}}$ denote the upper and lower temperature bounds of the storage system.

The remaining parameters were derived from system documentation and planning materials. The complete energy system model was used as the environment for the RL agent, which was implemented using Python and the Gymnasium framework [33], and pytorch [34]. The MILP model was formulated and solved using the Gurobipy interface to the Gurobi optimizer [35]. Both models rely on the parameters listed in Table 3.

2.4. Reinforcement Learning

This subsection outlines the deep reinforcement learning (RL) algorithm employed to the control of the heat pump and the heating element to reduce the energy costs and improve the self consumption. A Deep Q-Network (DQN) was implemented, following the approach presented in our earlier work [36]. The architecture of the DQN is illustrated in Figure 3.

The method builds upon the deep Q-learning algorithm originally introduced by Mnih et al. [37], and incorporates key improvements proposed by Van Hasselt et al. [38], which involve using two separate neural networks: a policy network for selecting actions and a target network for estimating Q-values. To enhance training stability, the algorithm integrates soft target updates as described by Lillicrap et al. [39]. Additionally, it uses experience replay with mini-batch sampling, following the stabilization technique as applied in the original DQN study [37].

In the present EC, the heating system’s controller acts as the agent, controlling the heat generators HP and HE, while the rest of the EC acts as the environment based on external conditions and actions. The agent’s control action consists of adjusting the SOC of the TES. The action set is discretized on a scale from 0 to 100, where 0 represents the minimum allowable TES charge and 100 the maximum, as regulated by the P-controller.

The optimization horizon covers one day, partitioned into 96 discrete time steps of ${\displaystyle \frac{1}{4}}$ h each. Consequently, one episode spans 96 steps. The system state at each step i is composed of the current ${\mathrm{SOC}}_i$ (scaled between 0 and 100), the number of remaining steps $(96-i)$ in the day, the forecasted electricity prices ${\pi }_j$ for the upcoming 96 intervals, the feed-in tariff f, the predicted thermal output ${\dot{Q}}_{\mathrm{out},j}$, the net electrical loads $P_{\mathrm{net}-\mathrm{load},j}$, and the seasonal indicator variables $B_1,B_2,B_3,B_4$ for those intervals. The forecasts for electricity prices and thermal load are assumed available one day in advance, reflecting realistic operational conditions where day-ahead market prices are published prior to execution. For the purposes of this application, operation is assumed to occur under perfect predictions without uncertainty, thereby simplifying the optimization problem. After each step, the forecast window shifts forward accordingly; unknown future values beyond the forecast horizon are set to zero. To enhance training stability, most state variables are scaled.

Electricity prices, used in calculating energy costs and rewards, were scaled according to a min-max normalization over the course of one day to a range from 1 to 10 via:

$${\pi }_i^{*}=9·{\displaystyle \frac{{\pi }_i-min\left(\pi \right)}{max\left(\pi \right)-min\left(\pi \right)}}+1$$

(13)

Similarly, the feed-in tariff f, which remains constant throughout each episode, was scaled using the same minimum and maximum values as the electricity prices to ensure consistent normalization:

$$f^{*}=9·{\displaystyle \frac{f-min\left(\pi \right)}{max\left(\pi \right)-min\left(\pi \right)}}+1$$

(14)

Scaling electricity prices serves three main purposes: first, it ensures that all state variables have comparable magnitudes; second, it normalizes daily price values to account for seasonal fluctuations; and third, it enables the agent to generalize to price scenarios that differ from those seen during training. Importantly, scaled prices were constrained to be greater than or equal to one, which avoids zero-cost intervals that might otherwise encourage unrealistic energy usage.

Thermal loads are scaled to represent the percentage of the TES capacity consumed per time step:

$${\dot{Q}}_{\mathrm{out},i}^{*}={\displaystyle \frac{{\dot{Q}}_{\mathrm{out},i}·\Delta t}{C_{\mathrm{stor}}(T_{\mathrm{ub}}-T_{\mathrm{lb}})}},$$

(15)

where $C_{\mathrm{stor}}$ is the thermal capacity of the TES.

To keep the state vector as small as possible, the electrical load $P_{\mathrm{load}}$ and PV generation $P_{\mathrm{PV}}$ have been summarized in the singular net load $P_{\mathrm{net}}$:

$$P_{\mathrm{net}}=P_{\mathrm{PV}}-P_{\mathrm{load}}$$

(16)

Collectively, the state vector at time step i is given by:

$$\begin{array}{cc}\hfill S_i=(& {\mathrm{SOC}}_i,96-i,{\pi }_i^{*},\dots ,{\pi }_{i+95}^{*},f^{*},{\dot{Q}}_{\mathrm{out},i}^{*},\dots ,{\dot{Q}}_{\mathrm{out},i+95}^{*},\hfill \\ \hfill & P_{\mathrm{net},i},\dots ,P_{\mathrm{net},i+95},B_1,B_2,B_3,B_4)\hfill \end{array}$$

(17)

The discrete action space from 0 to 100 corresponds to the target SOC (${\mathrm{SOC}}_{\mathrm{set}}$) for the P controller, from which the temperature setpoint $T_{\mathrm{set}}$ is computed as:

$$T_{\mathrm{set}}=T_{\mathrm{lb}}+{\displaystyle \frac{100-{\mathrm{SOC}}_{\mathrm{set}}}{100}}·(T_{\mathrm{ub}}-T_{\mathrm{lb}})$$

(18)

The environment dynamics were modeled using Equation (9) alongside a P-controller. Although each episode consists of 96 discrete time steps representing a full day, the P controller operates at a higher frequency within each simulation step, executing 60 control updates per time interval. The scaled price signal was then incorporated into the reward function, which distinguishes between grid consumption and feed-in:

$$R_i=\left\{\begin{array}{cc}-P_{\mathrm{grid},\mathrm{pos},i}·{\pi }_i^{*},\hfill & \mathrm{if}{P}_{\mathrm{grid},\mathrm{pos},i}>0\left(\mathrm{consumption}\right)\hfill \\ P_{\mathrm{grid},\mathrm{neg},i}·f,\hfill & \mathrm{if}{P}_{\mathrm{grid},\mathrm{neg},i}>0(\mathrm{feed}-\mathrm{in})\hfill \end{array}\right.$$

(19)

where $P_{\mathrm{grid},\mathrm{pos}}$ and $P_{\mathrm{grid},\mathrm{neg}}$ denote the electrical powers, consumed or fed in the grid, respectively, the ${\pi }_i^{*}$ denotes the scaled electricity price at step i and f the constant feed-in tariff.

The hyperparameters of the RL algorithm are summarized in Table A1. The following steps were taken for training the RL agent:

• At the start of each episode, a day index was randomly sampled from the training dataset.
• For the selected day, the corresponding normalized profiles were loaded, including electricity prices ${\pi }^{*}$, thermal load demand ${\dot{Q}}_{\mathrm{out}}^{*}$, electrical load $P_{\mathrm{net}}$, and feed-in tariffs $f^{*}$.
• Seasonal indicators such as winter ($B_1$), spring ($B_2$), summer ($B_3$), and autumn ($B_4$) flags were extracted for the selected day. This approach ensures diversity across episodes and captures a wide range of operational conditions.
• The agent was trained over a total of $5,000$ episodes.

2.5. Mixed Integer Linear Programming

In this section, the MILP model is introduced and described. The system was modeled and optimized to determine the global optimum of the EC’s operation. In order to ensure a fair comparison with the RL approach, the MILP model was provided with perfect prediction for each day of operation. This includes complete knowledge of the overall heat demand ${\dot{Q}}_{\mathrm{out}}$, electricity prices $\pi$, feed-in tariffs f, ground probe temperatures $T_{\mathrm{GP}}$, and photovoltaic generation $P_{\mathrm{PV}}$. The optimization was carried out over a horizon of $N=97$ discrete time stamps, resulting in 96 time periods p between them as $p\in \mathcal{P}=\{0,\dots ,N-1\}$, covering one full day.

The parameters used for the MILP formulation correspond to Table 3. The objective function minimizes the net electricity cost by considering the grid import power ($P_{\mathrm{grid},\mathrm{pos}}$), the electricity price signal $\pi$, the feed-in power ($P_{\mathrm{grid},\mathrm{neg}}$), and the feed-in tariff f, as well as the number of time steps N. The objective function is formulated as:

$$min\sum _{p=0}^{N-1}\Delta t\left(P_{\mathrm{grid},\mathrm{pos},p}·\pi {}_p-P_{\mathrm{grid},\mathrm{neg},p}·f_p\right)$$

(20)

This objective was calculated with a fixed horizon of 24 hours, more precisely defined as 96 quarter-hour intervals. Since each day was optimized independently and the total cost was calculated at the end, the following temperature bounds for $T_{\mathrm{stor}}$ were considered for consistency:

$$\begin{array}{cc}\hfill T_{\mathrm{stor}}\left(0\right)& =\left\{\begin{array}{cc}40\hfill & \mathrm{on}\mathrm{the}\mathrm{first}\mathrm{day}\hfill \\ T_{\mathrm{stor}}^{\mathrm{prev}}\left(N\right)\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill T_{\mathrm{stor}}\left(N\right)& =\left\{\begin{array}{cc}40\hfill & \mathrm{on}\mathrm{the}\mathrm{final}\mathrm{day}\hfill \\ \in [T_{\mathrm{lb}},T_{\mathrm{ub}}]\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(22)

As seen in the equations, the first day was initialized with a storage temperature $T_{\mathrm{stor}}$ of 40°C. On all following days, except the last one, the optimization began with the final temperature of the previous day, i.e., $T_{\mathrm{stor},N}^{\mathrm{prev}}$, to ensure that no energy is artificially lost or gained. While the optimizer was allowed to chose the end-of-day temperature of each day freely, on the last day the final storage temperature was fixed to $T_{\mathrm{stor},N}=40°\mathrm{C}$ to maintain consistency with the RL.

All other constraints of the model are listed below and must be satisfied at every timestep on each day. These include operational limits, system dynamics, and technical constraints, all of which ensure the physical and economic feasibility of the optimization results across the full time horizon.

$$\begin{array}{cc}\hfill \forall p\in \mathcal{P}:& a_p={\displaystyle \frac{hA}{C_{\mathrm{stor}}}}\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & b_p={\displaystyle \frac{{\dot{Q}}_{\mathrm{HE},p}+{\dot{Q}}_{\mathrm{HP},p}-{\dot{Q}}_{\mathrm{out},p}+hA·T_{\infty }}{C_{\mathrm{stor}}}}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & T_{\mathrm{stor},p+1}={\displaystyle \frac{b_p}{a_p}}+\left(T_{\mathrm{stor},p}-{\displaystyle \frac{b_p}{a_p}}\right)·e^{(-a_p·\Delta t)}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & {\dot{Q}}_{\mathrm{HP},p}=0.5·{\mathrm{COP}}_{\mathrm{Carnot}}·P_{\mathrm{HP},p}\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & {\dot{Q}}_{\mathrm{HE},p}=P_{\mathrm{HE},p}\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & P_{\mathrm{grid},p}=P_{\mathrm{grid},\mathrm{pos},p}-P_{\mathrm{grid},\mathrm{neg},p}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & B_{\mathrm{pos},p}+B_{\mathrm{neg},p}\le 1\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & P_{\mathrm{grid},\mathrm{pos},p}\le B_{\mathrm{pos},p}·M\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & P_{\mathrm{grid},\mathrm{neg},p}\le B_{\mathrm{neg},p}·M\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & P_{\mathrm{grid},p}+P_{\mathrm{PV},p}=P_{\mathrm{HP},p}+P_{\mathrm{load},p}+P_{\mathrm{HE},p}\hfill \end{array}$$

(32)

All variables used are bound by their respective boundaries as follows:

$$\begin{array}{cc}\hfill 1-1\forall p\in \mathcal{P}:& T_{\mathrm{lb}}\le T_{\mathrm{stor},p+1}\le T_{\mathrm{ub}}\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & 0\le P_{\mathrm{grid},\mathrm{pos},p}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & 0\le P_{\mathrm{grid},\mathrm{neg},p}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & 0\le P_{\mathrm{HP},p}\le P_{\mathrm{HP},max}\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & 0\le {\dot{Q}}_{\mathrm{HP},p}\le {\dot{Q}}_{\mathrm{HP},max}\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & 0\le P_{\mathrm{HE},p}\le P_{\mathrm{HE},max}\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill & 0\le {\dot{Q}}_{\mathrm{HE},p}\le {\dot{Q}}_{\mathrm{HE},max}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & B_{\mathrm{pos},p}\in \{0,1\}\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill & B_{\mathrm{neg},p}\in \{0,1\}\hfill \end{array}$$

(41)

As already mentioned in the description of the reinforcement learning approach, the 96 values forecasted for the electricity price $\pi$, feed-in tariffs f, overall heat demand ${\dot{Q}}_{\mathrm{out}}$, electrical load $P_{\mathrm{load}}$, and solar generation $P_{\mathrm{PV}}$ were assumed to be perfectly accurate under the assumption of perfect prediction. This assumption simplifies the optimization by neglecting forecasting errors, allowing the model to focus solely on the system’s operational optimization. However, in practical applications, forecasting uncertainties do impact the solution quality and must be considered for a more robust approach.

3. Results and Discussion

In this section, the results of both optimizers are presented and discussed. In Figure 4, various electric energy flows P within the system are presented. Panel a displays five days of electrical power flows under the RL-optimized operation, while panel c shows the behavior of the MILP-based system. Power flows into the hub, including photovoltaic generation ($P_{\mathrm{PV}}$) and grid import ($P_{\mathrm{grid},\mathrm{pos}}$), are represented as positive values. In contrast, power flows leaving the hub, such as the electrical load ($P_{\mathrm{load}}$), the heat pump ($P_{\mathrm{HP}}$), the heating element ($P_{\mathrm{HE}}$), and electricity exported to the grid as feed-in ($P_{\mathrm{grid},\mathrm{neg}}$), are represented as negative values. The corresponding real-time electricity price signal and the constant feed-in tariff are shown in panel b.

Although the general operation patterns appear similar at first glance, several key differences emerge between the two control strategies, most notably in the operation of the auxiliary heating element. For the heating element, the RL agent tends to favor longer periods of moderate operation or to avoid its activation altogether. In contrast, the MILP model opts for short, high-intensity bursts, operating the auxiliary heater at full capacity when required. This indicates a more aggressive but shorter-duration heating strategy under the MILP control to achieve the global optimum of operation. Both optimization approaches clearly attempt to maximize the self-consumption of PV-generated electricity before feeding into the grid. The constant feed-in tariff of 0.08 €/kWh appears insufficient to incentivize export when local usage is possible. Furthermore, the five-day excerpt illustrates a period of high thermal demand. Both systems frequently operate the heat pump and auxiliary heating elements at or near their maximum power levels, indicating a substantial heating demand driven by cold weather conditions (below -4°C) for the EC. These conditions correspond to the second percentile of all annually occurring ambient temperatures.

Figure 5 depicts the transient energy balance of the TES depicted as a time-series for the same time period as depicted in Figure 4. It showcases the RL-controlled and the MILP-controlled operations, panel a and panel c, respectively. Positive flows symbolize heat entering the sensible storage tank during the operation of the heating element or the heat pump; negative flows therefore depict the heat taken from the TES, i.e., the heating load ${\dot{Q}}_{\mathrm{out}}$ for space heating. The sum of the heat flows is depicted as ${\dot{Q}}_{\mathrm{charge}}$ and reflects the charging power (or discharging power, if negative) of the TES within a period. In panel b, the temperatures of both operation modes are shown in degrees centigrade as well as in % SOC.

Starting with the temperature curves in b, the general shapes are quite similar, indicating a similar approach by both RL and MILP. The main difference between the two is that RL manages to dip under the lower temperature constraint, which prompts the logic of the RL to activate the heating element. As the MILP is bound by the temperature constraints, it is not allowed to let this happen. For the general operation seen in a and b, the approach seems quite similar, as both systems have the same heat load applied to them. For these five days, the RL seems to have a more bipolar behavior, not using the inverter function of the heat pump, but instead turning it on and off, while the MILP operation seems to utilize the modulation ability of the heat pump, resulting a more continuous and moderated operation of the storage tank, while the RL-optimized TES exhibits more frequent and abrupt switching between charging and discharging.

To compare the performance of the reference scenario with the RL and MILP optimization approaches, five key performance indicators (KPIs) were selected, as illustrated in Figure 6. Panel a presents the relative deviation of both optimization methods from the reference case in terms of cumulative costs, ${\mathrm{CO}}_2$eq emissions, grid power consumption, self-consumption ratio (SCR), and self-sufficiency ratio (SSR). In addition to these relative deviations, Table 4 provides the corresponding absolute values. Panel b displays the cumulative cost evolution over the 105-day testing period, highlighting the comparative behavior of the MILP and RL approaches.

As seen in panel a, MILP and RL are generally quite similar in their results. Both optimizers achieve cost savings compared to the reference scenario. While the MILP is able to find the global optimum at 10.06% reduction compared to the reference scenario, the RL optimization is able to approach this with 8.75%. Although not optimizing for emissions, ${\mathrm{CO}}_2$ emissions are also reduced in the RL and MILP case, by 4.94% and 4.01%, respectively. Grid power usage is nearly the same irrespective of the optimization strategy. One reason why emissions could be reduced is the increase in the SCRs and SSRs. As solar power was considered as ${\mathrm{CO}}_2$-neutral, an increase in self-consumption leads to a reduction in ${\mathrm{CO}}_2$. This also explains the better results of the RL optimization, as it achieved the highest self-consumption rate, thereby also reducing ${\mathrm{CO}}_2$ emissions.

The cumulative cost curves in panel b reveal that the overall shape remains similar across both optimization approaches and the reference scenario. However, starting in spring, the curves begin to diverge as the optimized strategies achieve increasing cost savings. This divergence is driven by two opposing trends: higher PV production yields more energy that can be used directly, while milder temperatures reduce the heating load and thus the overall energy demand . Over the full year, MILP and RL achieve total savings of 67.58 € and 58.76 €, respectively, with MILP delivering the best overall result and RL performing comparably close.

4. Conclusions

This study shows that optimization algorithms like RL can approach the global optimum as determined via MILP to a degree where the difference is negligible in practice. Based on the measured KPIs, RL even outperforms MILP in some areas, including CO₂, SCR, and SSR, making it a strong candidate for real-world applications. This study distinguishes itself by using high-resolution, real-world data, including measured demand, PV generation, and real-time electricity prices, under near-identical conditions for both methods.

Even if the result is promising and the optimization itself is faster in RL than it is for MILP, there are some aspects to consider. RL has two main downsides. First, a large amount of data is needed to train and validate the agent, which even required the implementation of additional features such as seasonality or normalized data. Secondly, training for $5,000$ episodes took multiple hours to complete, meaning that from point zero to optimization, MILP was much faster than RL, the reason being the multidimensional state of each training episode combined with the relatively fine temporal resolution of 15-minute steps. Another topic which cannot be discounted is the implementation of the auxiliary heating element in the RL. The kind of implementation used allows the RL-optimized system to drop below the hard lower bounds of the MILP by two degrees before applying the heating element, leading to some temperature violations, which gives RL slightly more flexibility overall. However, a closer inspection shows that these violations are generally mild, with only 1 % of cases falling below 34.0°C and none below 33.0°C.

Despite these limitations, RL demonstrates strong performance under realistic conditions. With some targeted preprocessing, it is able to achieve near-optimal results and offers the key advantage of real-time control once training is complete. In deployment, the agent can act autonomously and adapt over time through retraining or transfer learning, making RL a promising solution for dynamic environments and potentially suitable for edge applications. Both approaches should be able to improve the performance of a simple system, so the question might not always be which method can produce the better result, but rather which method is more feasible for real-world application, as small and fast computers become increasingly available.

Future research could focus on validating these findings in real-world settings, for example through lab setups using synthetic loads. Expanding the system to include additional components like solar thermal, biomass combined heat and power, or decentralized storage would test both scalability and robustness. Another valuable direction is the integration of real load forecasting and resilience to unexpected events such as sudden demand shifts or system faults. Beyond the immediate results, this work contributes to a growing body of evidence that learning-based control strategies can approach traditional optimization methods under certain real-world constraints. The insights gained from this study can be applied to similar energy systems and extended to other flexibility assets, supporting the transition from academic methods to real-world deployment.

Ultimately, this work delivers one of the first direct, high-resolution comparisons of RL and MILP for optimizing an existing EC using real-world data, PV generation, and market signals. It demonstrates that deep RL is not only a viable alternative to MILP but can in some cases even match a carefully tuned MILP model, while offering practical advantages in speed, adaptability, and scalability for real-time control of complex energy systems.