Two-stage stochastic, large-scale optimization of a decentralized energy system: a case study focusing on solar PV, heat pumps and storage in a residential quarter

Abstract

The expansion of fluctuating renewable energy sources leads to an increasing impact of weather-related uncertainties on future decentralized energy systems. Stochastic modeling techniques enable an adequate consideration of the uncertainties and provide support for both investment and operating decisions in such systems. In this paper, we consider a residential quarter using photovoltaic systems in combination with multistage air-water heat pumps and heat storage units for space heating and domestic hot water. We model the investment and operating problem of the quarter’s energy system as two-stage stochastic mixed-integer linear program and optimize the thermal storage units. In order to keep the resulting stochastic, large-scale program computationally feasible, the problem is decomposed in combination with a derivative-free optimization. The subproblems are solved in parallel on high-performance computing systems. Our approach is integrated in that it comprises three subsystems: generation of consistent ensembles of the required input data by a Markov process, transformation into sets of energy demand and supply profiles and the actual stochastic optimization. An analysis of the scalability and comparison with a state-of-the-art dual-decomposition method using Lagrange relaxation and a conic bundle algorithm shows a good performance of our approach for the considered problem type. A comparison of the effective gain of modeling the quarter as stochastic program with the resulting computational expenses justifies the approach. Moreover, our results show that heat storage units in such systems are generally larger when uncertainties are considered, i.e., stochastic optimization can help to avoid insufficient setup decisions. Furthermore, we find that the storage is more profitable for domestic hot water than for space heating.

Appendices

Appendix A: Further information on the energy supply and demand profiles

Concerning the supply side profiles, we provide details on the validation of the global solar radiation output of the developed Markov model (see Section 3.1) for illustrative purposes. Moreover, we present information about how the radiation profiles differ between the different scenarios. The other Markov model output parameters (temperature and cloudiness) have been validated accordingly. Note that historical measurement data over a longer horizon are only available in hourly resolution. For the validation, we therefore aggregate the model output from 15 min to hourly resolution. Let $\rho _a =\left( {\rho _a^1 ,\,\ldots ,\rho _a^{8760} } \right)$ be an hourly series of global solar radiation in year a, where $a\in \left\{ {1971,\,\ldots ,\,2011} \right\}$ for the historical data and $a\in \left\{ {1,\,\ldots ,\,100} \right\}$ for the results of the Markov model. We now validate the Markov model on the basis of four indicators: (1) the total annual radiation supply in year a defined by $\hbox {P}_a =\mathop {\sum }\nolimits _{h=1}^{8760} \rho _a^h$ as a long-term indicator and (2) the hourly volatility in year a defined by $\mathrm{vola}_a =\sigma \left( {\rho _a } \right) /\mu \left( {\rho _a } \right)$ as a short-term indicator, where $\sigma \left( {\rho _a } \right)$ and $\mu \left( {\rho _a } \right)$ are the standard deviation and arithmetic mean of the global solar radiation in year a, respectively. In addition, we consider (3) the maximum amplitude of radiation supply $(\mathrm{MARS}_{a})$ in year a and (4) the maximum gradient of radiation supply $(\mathrm{MGRS}_a )$ in year a as defined by Schermeyer et al. (2015). In order to validate the performance of the Markov model in the long run, we compare the arithmetic means $\mu \left( \cdot \right)$ of these four indicators over all available years (simulation results vs. historical data) as well as the 5 and $95\%$ quantiles $\mathrm{quant}_{5\% } \left( \cdot \right)$ and $\mathrm{quant}_{95\% } \left( \cdot \right)$ over all available years to analyze the range of variation. Table 3 shows the relative deviation of these means and quantiles of the four indicators between the Markov model results and historical data. For instance, the values in the column $\mu \left( \cdot \right)$ are calculated as $\left( {\mu ^\mathrm{Mod}\left( \cdot \right) -\mu ^\mathrm{Hist}\left( \cdot \right) } \right) /\mu ^\mathrm{Hist}\left( \cdot \right)$, where the superscript Mod denotes the model results and the superscript Hist denotes historical data. The values in the other columns are calculated accordingly. Overall, this comparison shows satisfying results. At the same time, however, the table shows that there is room for further improvement of the Markov model in the future.

Table 3 Results of a comparison of simulation model results and historical weather data for four indicators

Figure 8 shows the importance of using 15-min profiles rather than hourly profiles on the supply side. The left diagram shows the variability of PV power output between the 100 considered scenarios in general for a day in June. It also shows that the spikes on the top only occur during short periods of time (15-min intervals rather than hours). This implies that the maximum amplitude of radiation supply is underestimated with hourly profiles. Moreover, when it comes to choosing the optimal sizes of the energy system’s components, the gradients of power output between time steps are very important. This is particularly relevant for storages which are at the core of our case study. The right diagram in Fig. 8 shows that the maximum positive and negative gradients are strongly underestimated when using an hourly resolution as opposed to a 15-min resolution.

Figure 9 shows the optimal storage size for SH and DHW under SMILP-1 and SMILP-2 for different temporal resolutions. When time steps of 60 min are used instead of 15 min, the optimal storage sizes differ by up to $50\%$. The lower temporal resolution reduces the load-shifting potential and leads to smaller storage units for SMILP-1 (without integers at the second stage). In case of SMILP-2, the stepwise flexibility of the heat pumps is reduced when moving from a 15-min resolution to 60 or even 120 min. This makes the storage units more attractive. This effect outbalances the reduced load shift potential and results in larger storages for SMILP-2 for coarser temporal resolutions.

Fig. 8

Variability of PV generation between scenarios (left) and range of maximum gradients (right) for an illustrative June day

Fig. 9

Influence of the temporal resolution on the optimal storage size for space heating (SH) and domestic hot water (DHW)

Concerning the demand-side profiles, as described in Sect. 3.2, we use a standard load profile approach (based on so-called H0 profiles), to generate electricity demand time series for the 70 households of the quarter. Thereby, the total yearly electricity consumption (without the electricity demand of the heat pumps) is calculated according to VDI 4655, which takes the number of residents and the usable floor surface into account. Aiming for an ex post validation of the assumption that 70 households can be approximated by H0 profiles, we compare the H0 profiles to measured electricity demand profiles of households that have already moved into their dwellings in the quarter (see Table 4). The comparison is based on 40 households since we only include households where measurement data are available for an entire year and the remaining households have moved in at a later date. Table 4 shows the (linear) correlation coefficient between the H0 profiles and the measured profiles, the mean absolute percentage error (MAPE), the root-mean-square percentage error (RMSPE) and the relative difference of the demand volatility $\Delta \left( \mathrm{vola_a } \right)$, where $\Delta \left( \mathrm{vola_a } \right) =\left( \mathrm{vola_a \left( {H0} \right) -\mathrm{vola}_a \left( \mathrm{measured\,profiles} \right) } \right) /\mathrm{vola}_a \left( \mathrm{measured\,profiles} \right)$. The correlation coefficient between the 40 households and the H0 profiles already amounts to $78\%$. For larger numbers of households, Hayn et al. (2018) show that the correlation coefficient between 100 households and the H0 profile increases to $90\%$. We therefore expect the correlation coefficient of the entire quarter to be between 78 and $90\%$. In terms of the load volatility, we find that there is only a $-7\%$ difference between our 40 households and the H0 profiles, which we expect to further decrease for 70 households (similar to the effect described for the correlation).

Table 4 H0 profiles compared to real-measured electricity demand of 40 households of the quarter

Table 5 Possible DFO algorithms for the outer optimization of the fixed variables (based on Rios and Sahinidis 2013; Conn et al. Conn et al. 2009)

Appendix B: Overview of derivative-free optimization approaches

As described in Sect. 3.3, a derivative-free optimization (DFO) is used for the outer optimization of the fixed mixed-integer variables. Table 5 summarizes DFO methods with regard to the mentioned requirements. DFO refers to problems when information on the derivatives of f is unavailable, unreliable or impractical to obtain. This definition includes any algorithm applied to these problems, even if the algorithm involves the computation of derivatives for functions other than f (Rios and Sahinidis 2013). Not included are commonly known algorithms such as branch-and-bound, cutting plane or Lagrangian relaxation, albeit specific variants could be considered as DFO according to this definition. However, the economic optimization subsystem introduced in Sect. 3.3, particularly the approach to parallelization, would need to be substantially changed to apply such algorithms.

The main task is to determine optimal storage sizes of the residential quarter. Therefore, the problem is decomposed by fixing storage sizes and storage levels. Since the non-optimization of the storage levels leads to a negligible error, only few fixed integer variables need to be optimized, i.e., 8 storage sizes of the quarter. To this purpose, SAHC is sufficient for the outer optimization. The advantages of the implementation are its simplicity, flexibility and reliability. Furthermore, it robustly proceeds to the (local) optimum even with inaccurate solutions of the subproblems. The complete SAHC procedure is presented in the following:

Procedure of the steepest-ascent hill-climbing (SAHC) method:
Step 0:	(Initialization) compute $f\left( \varphi \right)$ for an initial $\varphi$ (e.g., $\varphi =0$) by using POP and set step size $s_k$ for each fixed variable $\varphi _k$ of vector $\varphi$. If $\varphi _k \in \mathbb {Z}$, then $s_k \in \mathbb {Z}$. Let $e_k \in \mathbb {R}_+^{I+R-v}$ be the k-th unit vector, where I is the number of fixed first-stage variables and $R-v$ is the number of fixed second-stage variables.
Step 1:	Add $s_k$ to $\varphi _k$ and compute $f\left( {\varphi +s_k e_k } \right)$ and subtract $s_k$ from $\varphi _k$ and compute $f\left( {\varphi -s_k e_k } \right)$ by using POP sequentially for each fixation $1\le k\le I+R-v$. Note if $f\left( {\varphi +s_k e_k } \right) >f\left( {\varphi -s_k e_k } \right)$, then $step_k^*=+s_k e_k ,$ else $step_k^*=-s_k e_k$.
Step 2:	Select $\varphi ^{*}\in \left\{ {\varphi \pm s_k e_k |\forall 1\le k\le I+R-v} \right\}$ with $f\left( \varphi \right) =\mathop {\hbox {min}}\limits _k \left\{ {f\left( {\varphi \pm s_k e_k } \right) } \right\}$.
Step 3:	Define $\Delta f\left( \varphi \right) _{rel} =\left( {f\left( \varphi \right) -f\left( {\varphi ^{*}} \right) } \right) /f\left( \varphi \right)$.
Step 4:	If $\Delta f\left( \varphi \right) _{rel} \quad \le \quad 0$, then $s_k =\frac{s_k }{2}$; if $\varphi _k \in \mathbb {Z}$ and $\frac{s_k }{2}<1,$ then go to step 6; if $\varphi _k \in \mathbb {Z}$ and $\frac{s_k }{2}\notin \mathbb {Z}$, then round $\frac{s_k }{2}$ to the larger integer; go to step 1. Otherwise continue.
Step 5:	If $\Delta f\left( \varphi \right) _{rel}>$ stopping criterion $a\in \mathbb {R}_+$, then accept $f\left( \varphi \right) =f\left( {\varphi ^{*}} \right)$ and $\varphi =\varphi ^{*}$, compute $f\left( {\varphi +step_k^*} \right)$ by using POP sequentially for each fixation $1\le k\le I+R-v$ and go to step 2. Otherwise continue.
Step 6:	(End) Stop. The local optimal solution value is $f\left( {\varphi ^{*}} \right)$ with the vector $\varphi ^{*}$.

Appendix C: Further information about the mathematical model of the quarter

Table 6 lists the complete nomenclature of the residential quarter modeled as a two-stage stochastic mixed-integer program.

Table 6 Nomenclature