An overview of the whole approach is presented in Figure 1. The main novelty of HPSS is the use of the subscripting feature, i.e., an array structure of elements that can be repeated as many times as the modeler needs. Taking advantage of this feature, the dispatch of electricity is configured through variables containing subscripts of 8760 elements (the hours in a year). These are used to import the hourly profiles of demands and VRES, calculate the intermediary steps, and report hourly indicators of interest, such as curtailment^† or power ramps. Concurrently, the integration method solves the differential equations that represent the dynamics of the capacities of technologies, the evolution of efficiencies, and the aggregated indicators over time. Consequently, subscripts have generated the ability to represent two or more discrete temporal dimensions while the continuous integration process is ongoing. This idea is not new and has been employed in other dissimilar fields facing similar problems, e.g., irrigation [25] and migration [26].

Figure 1.

General overview of the HPSS model. Yellow variables are calibrated for the base years. Grey boxes are exogenous data; while blue boxes are exogenous data in this work, but endogenous in an IAM (to be endogenized in further work). The equations set out in this work are summarized as “Addition” or “Subtraction” in the gray color. The two intermittent boxes are variables repeated to clarify the diagram. It is well known that “the system must be able to meet the continually changing load demand” [28], verifying the first property of the allocation function.

Nevertheless, the question about the dispatch of electricity that needs to be answered is: How can the load demand be fulfilled from different suppliers? Here, the one-to-many allocation function takes place. It is a mathematical formulation that solves a single market, where several suppliers are trying to satisfy the demand. In this work, the functionality reproduces the operation of the power system every hour. The software used to develop the model is Vensim DSS, where the “ALLOCATE AVAILABLE” function is included in the software set [27]. In short, the allocation function has six properties implicitly assumed and specified for our problem:

• Conservation of matter: the amount received by the demanders has to be equal to the amount managed by the suppliers.

• Nonnegative values: all the quantities must be equal to or greater than zero.

• Conservation of intent: no supplier shall provide more than the desired/possible amount of energy to supply. Similarly, no demander shall receive more than what has been chosen.

• No loopholes: under adequate conditions, each agent should receive/supply its stated demand/offer.

• Differentiation: in the case of energy shortage, extremely low priority agents should receive/supply little or nothing and extremely high priority agents should receive/supply a higher amount of energy or almost everything.

• Continuity: with small changes in priorities, supply and demand should cause small and smooth changes in the resulting allocations.

The non-dispatchable generation (VRES), nuclear (N), imports (I), and exports (E) load is subtracted from the load demand (D) in order to place the difference (residual load) in either a positive or negative distribution (RL). This concept of residual load (RL) [29] has been used to fulfill the second property throughout equation (1) (where h corresponds to each hour). Positive values (PRL) correspond to hours with a deficit of energy before management; while negative values mean over-productive hours (an excess of electricity). Both resulting positive and negative distributions are treated as separate markets to define which technologies can respond by covering the deficit (e.g., thermal power plants), increasing the load demand (e.g., heat pumps), reducing the generation (e.g., curtailment), or influencing both markets (e.g., storage).

$R{L}_{\text{h}}=D_{\text{h}}-VRE{S}_{\text{h}}-N_{\text{h}}-I_{\text{h}}+E_{\text{h}}\text{ 1}\le \text{h}\le \text{8760}$ (1)

The maximum load of the suppliers is limited by the installed capacity of the technology in the year (IC) and, in the case of hydropower units, the maximum hourly producible^‡ profile.

On the other hand, every dispatchable power plant has a range of technical flexibility to offer an hourly potential generation (G_h) due to the inertia of the heating system and the turbines, thermal stress and security risks, among others (section 2 in [30]). So the generation of electricity is partially restricted so as to reproduce the adaptability of the system following the demand, an effect also related to the schedule that suppliers create with the objective of maximizing profits. This effect is reproduced with equation (2) in the model, which divides the generation into a flexible load and an invariable load controlled by a flexibility parameter (fp). An exception has been introduced for the case of storage, which is explained later. This method respects the third and four properties of the allocate function that is applied later.

$G_{\text{h}}=IC\times \left(1-fp\right)+IC\times fp\times \frac{PR{L}_{\text{h}}}{\underset{1\le h\le 8760}{{\displaystyle \max }}\left(PR{L}_{\text{h}}\right)}$ (2)

Competition for the demand load is present in any liberalized electricity market. A priority according to the technology is assumed in order to manage the energy shortages in hours when the aggregated desired supply is higher than the load demand. This effect is property number 5 of the allocate function. In this work, the priority is an indicator of how good the position of a supply technology is in the electricity market. The higher the priority, the higher the energy delivered. Since the suppliers cannot deliver more energy than the demand, a smaller load implies more competition. Finally, the model assumes that small modifications in priorities deliver small changes in the resulting supply mix; given that the hour, as a time unit, is sufficiently short as to assume smoothness in the energy delivery. This continuity is achieved by applying an exponential function to the shape of the priority response for the supply curves (green line in the graph of Figure 1).

Pumped hydropower energy storage (PHES) has been modeled in a singular way to avoid simultaneous equations in the code. From the literature, the knowledge states that: a) these units pump water during off-peak hours (low demand) and turbine water afterwards during peak-hours [31]; b) they cover from daily to seasonal cycles [31]; c) they can offer primary and secondary ancillary services to obtain profits, added to those obtained from the intraday and day-ahead markets (increasing profits compared to traditional bidding [32], even double [33]). The behavior of PHES becomes even harder in potential smart energy systems prepared for a high, or even 100%, renewable share in the energy system, where more technologies may offer similar services (vehicle-to-grid) and competition for the renewable overproduction [34]. Representing the totality of this complexity falls beyond the scope of this work; however, some operational details have already been integrated to show the potentialities and barriers in modeling.

First, an estimation of the generation mix, which we called hourly pre-dispatch, is calculated by an allocation function to get the share of the positive residual load pumped to the storage from the market competition (PHES_pos_turb). As a consequence, the annual consumption of the pumping mode is calculated with the efficiency of this technology (eq. (3)). Then, this consumption is split into 8760 hourly loads based on its historical profile (PHES_pos_cons).

$PHES\_pos\_con{s}_{\text{h}}=\frac{PHES\_pos\_tur{b}_{\text{h}}}{PHES\_efficiency}$ (3)

The consumption of the pumping mode is added to the load demand for a second allocation (dispatch) where, since it has already been accounted for, the turbine mode of PHES is not involved (LD_2nd_al). In short, the pre-dispatch is used as a kind of forecast to depict the turbine and pumping patterns of PHES. This pre-calculation avoids simultaneous equations.

As for the negative residuals (excess of electricity), storage is the only technology able to make a profit in this version of the model, a clear advantage in scenarios with high shares of non-dispatchable contribution. Thus, PHES acts as a monopolistic agent, pumping as much energy as possible (PHES_neg_turb). In this case, the consumption during the pumping mode does not increase the base load demand in hours with a deficit of energy, since the negative value is potentially curtailment. On the contrary, this storage technology transfers energy from the negative to the positive distribution, so the positive value of the load demand is reduced by a quantity equal to the overproduction gathered by storage (PRL2_h in eq. (4)).

$PRL{2}_{\text{h}}=PR{L}_{\text{h}}+PHES\_pos\_con{s}_{\text{h}}-PHES\_pos\_tur{b}_{\text{h}}-PHES\_neg\_tur{b}_{\text{h}}$ (4)

Therefore, the remaining negative residual loads (after storage) are accounted for as curtailment. Although there are other options to flexibilize the supply and demand sides facing deeper decarbonization marks, only storage has been considered due to the lack of data. Examples of other flexibility options are power-to-heat, power-to-hydrogen, or smart electric vehicles (section IX-F and throughout the text, respectively, in [34]). However, these options are residual in the historical period and data are not fully available (e.g., hourly profiles for calibration). The inclusion of more technologies in the competition for the excess of electricity is further work for the approach.

Once storage has been calculated, the second allocation re-calculates the electricity mix for the dispatchable technologies (LD_2nd_al). The priorities are the same for both allocate functions; they do not change from the pre-dispatch to the dispatch allocation or across hours.

Once the properties have been contextualized for the problem, the following paragraphs set out the technological specifications for the modeling of the Spanish power system. The dispatchable power plants (combined heat and power, CHP, and thermal power plants, PP) are sorted, in turn, by fuel: gas, geothermal, liquid, and solid fuels; as well as two hydropower technologies: traditional hydropower plants (dammed and mixed) and pump hydroelectric energy storage. In total, there are 12 technologies operating when necessary to fulfill the positive residual loads (Table 2), with a degree of freedom in the allocation mechanism, i.e., the hours when non-dispatchable generation is not enough. Table 2 also shows the coherency between databases and the HPSS model. The priorities are restricted to between 0 and 10, that is, the annual competitiveness in accessing the electricity market.

Hourly generation profiles of VRES, nuclear and international transactions are exogenous inputs introduced due to the high uncertainty around them. Forecasting techniques have been developed for different time frames. However, the behavior of solar irradiation and wind speed predictors remains unclear beyond 2 days ahead [35], while climate change also impacts the availability of resources in the long term [36][37]. Finally, international power exchanges and nuclear power plants are running under strict bilateral agreements to exchange electricity between different markets, making their profiles harder and unclear. The dynamics of human-decision-based generation, as well as the dynamization of the VRES profiles, are beyond the scope of this work. As mentioned later, the exogenous profiles introduced to calibrate the historical period are also used as input to select the patterns for the future context in which energy policies are to be applied.

A common way to estimate the generation is from the installed capacity. This is done via normalization. For the exogenous technologies, the hourly load is divided by the installed capacity. So, the capacity expansion proportionally determines the generation in the scenario. In equation (5), the correction parameters (CP) translate these installed capacities into hourly generations. When these values remain very close to one, it means that the databases are consistent with the hourly profiles introduced. CP values are estimated via calibration.

$G_{\text{year}}={\displaystyle \sum _1^{8760}{G}_{\text{h}}\times CP\times IC}$ (5)

The independence across elements of the subscript is a limitation of the approach, i.e., an element cannot be updated or calculated based on another element. Consequently, the energy mix of one hour cannot be used to calculate the following one. This implies that the turbine/pumping mode of hydroelectric plants, for instance, cannot be constrained by the level of water in the up and down reservoirs.

However, this limitation is partially solved after the allocate functions. The maximum power ramps in the year are estimated from the hourly generations so as to deliver an indicator of the necessary technical requirements, which may validate the scenario or not. A similar reasoning may be followed for the hourly exchanges of water.

The electricity demand follows the Eurostat’s convention on energy balances. The data are exogenously imported from 24 economic sectors, including the energy sector; which are then sorted into the extraction of fossil fuels (coal, natural gas, oil), consumption to produce coke oven coke for the steel and iron industry, and petroleum refineries. The transmission and distribution losses increase the electricity demand up to the demand at the busbars of the generators. A parameter multiplies this last demand to meet the demand from all databases (ParCalDatabase). As is shown in the next subsection, this multiplier matches the annual demand provided by the energy balances of Eurostat and the value reported by the power system operator, REE.

This section sets out the data sources used to study the Spanish power system with HPSS.

The annual consumption of electricity by the economic sectors and transmission losses have been gathered from Eurostat [38]. Moreover, the hourly profiles of the technologies have been collected from the power system operator’s database (ESIOS, [39]). Together, they make up the main data source.

Water resources, evaluated as producible energy, are extracted from the AEMET reports on a daily basis [40]. The estimation of the hourly potential for dammed and mixed hydropower generation is homogenously distributed over the day, so the daily value is divided by 24 hours.

Additionally, the data related to combined heat and power (CHP) generation has been found in a complete report written by the Institute for Energy Diversification and Saving (IDAE, Spanish acronym) [41]. Facilities are categorized by technology, sector, and fuel, including the installed capacity, electricity production, fuel consumption, and net heat production. There is a lack of information about residual generation in Spain from marine, solid biofuel, renewable municipal waste, and liquid biofuel.

Finally, the power system operator (REE) annually delivers reports [42] with Excel files that have been used for calibration purposes. For example, electricity generation by technology, including the net international balance and pumping hydropower consumption.

The availability of data in Spain has allowed us to study the calibration of most of the parameters over four years, from 2017 to 2020. This period is characterized by dissimilar meteorological conditions, the deployment of solar and wind (38%), the decommissioning of coal power plants (40%), and the disruptive effect of the Covid-19 disease in 2020.

Unfortunately, geothermal, run-of-river, oceanic, solar urban and offshore wind generation are negligible in Spain, there are no data about their hourly profiles. Consequently, they were not considered.

Powell’s algorithm is a hill-climbing method (iterative) that finds a function’s minimum value by sequential one-dimensional searches using the concept of conjugate directions [43]; the lack of gradients in this method is useful for non-differentiable and differentiable functions [44]. The Powell method was selected to calibrate and manage the priorities of the power system. This method is a Vensim software tool, explained in the documentation [27] and in specific publications [45][46]. The optimal values of priorities that minimize the payoff function are specified in Table 3; while the results of the calibration are shown in Table A. 1. As mentioned above, the priorities are constrained between 0 and 10 for the calibration of the allocation functions.

The parameter that evaluates the consistency between Eurostat and the power system operator database is very close to one in all the years (ParCalDatabase). So, there is perfect consistency between official agencies.

The hypothesis of assuming equivalence between the hourly average power and hourly generation is also acceptable, given the closeness to one in the results of international trade, wind, solar, and nuclear generation. Therefore, the estimation of the yearly generation from the hourly power profile and the installed capacity of the technology is also accepted.

However, the dimensionality of the problem has been simplified, i.e., the number of degrees of freedom in the calibration algorithm, in such a way that the maximum potential load generated by the dispatchable power plants corresponds to the installed capacity (fp = 0). Further analysis is required to explain the constraints in the generation of these technologies.

Figure 2 illustrates the priorities of dispatchable power units. Dam hydropower and CHP fueled by gas units are at the top of the electricity market in Spain, while the pumping mode of hydropower and PP units (traditionally run by fossil fuels) are at the bottom. The range of dots is narrower over time, which indicates that the market evolved towards a more competitive scenario. In the ranking, there is a significant change between 2018 and 2019, when PP solids (mostly coal) experienced a depletion in the installed capacity, promoting the role of other technologies, such as PP gas, CHP liquid and PHES, to replace it. Finally, the disruption of the Covid-19 disease did not affect priorities, indicating that these are driven by technical reasons of the suppliers rather than the electricity demand pattern.

Figure 2.

Priorities of dispatchable power plants in the historical period

The historical round-trip efficiency in PHES units was found to be between 54.3-61.2%, far below the common range of values in the literature (70-80%, [31]). As for the power network, the technical and non-technical transmission and distribution losses^§ have increased over the historical period from 9.8% to 11.4% (Eurostat). The reasons behind this increasing inefficiency are beyond the scope of this article.

Equation (6) accounts for the error in all the results of interest. Percentages between the actual and the calibrated annual generations are reported in Table A. 2, where all of them fell far below 5%. “CHP waste” delivered 4.04% more energy than expected in 2020; however, its generation is very low in relation to the total amount supplied (around 0.07%). The deviations in the share of renewables are +0.33% (2017), -0.09% (2018), -0.37% (2019), and -0.96% (2020). The highest error for PHES was only -1.21 MWh (2019). In agreement with the Spanish power system operator, PHES was considered a non-renewable energy source. Furthermore, when a set of technologies is mixed in the same category of the model, the installed capacity of each technology is used as a proxy to make the generation share, e.g., biomass represented 15.57% of the “PP solid fuels”, while biogas 1.01% of “PP gas fuels”, in 2020.

$Annualerror\left(\%\right)=\frac{realproduction-calibratedproduction}{realproduction}\times 100$ (6)

Table A. 3 shows acceptable errors for the purpose of the model. As can be seen, most of the means of the errors are very close to zero, apart from PP_solid & PP_waste, probably caused by the lack of disaggregation of these technologies in ESIOS. In short, the low errors in the model parameters validates the accuracy, given the available knowledge and data.

The rest of this subsection tests the behavior of the model when facing changes in two sensitive variables: First, the installed capacities of nuclear power plants; and second, the priority of PHES. In order to better visualize non-linearities in the results, the expansion of installed capacities and priorities follows a future linear trend. For the same reason, the demand remains constant over time and equal to the corresponding value of the historical year (268 TWh, 269 TWh, 265 TWh, 250 TWh, respectively).

The first exercise phases out the nuclear power plants by 2030. Meanwhile, the priorities remain constant and equal to the corresponding historical year. Four future trends are obtained, shown in Figure 3. The brown vertical line separates the historical period from the future.

Figure 3.

Evolution of some suppliers for different historical conditions. One color by year and one type of line by technology. The brown vertical line separates the historical period (left, one known point per year) from the scenario (right, four possible points per year)

The profile for nuclear power plants is exogenously assumed, and is therefore inflexible. This situation creates a reaction of dispatchable power plants, increasing their generation to fulfill the gap. The year 2017 was a dry year; therefore, hydropower generation barely grew over time (dashed blue line), being substituted by fossil fuels in thermal power plants (blue dash-dot line).

The political-medical decision for a lockdown due to the Covid-19 disease (2020) provoked a drastic reduction in the electricity demand, especially in households. It implied less energy consumption from fossil fuel power plants, but also higher hydropower generation, given the favorable hydrology of that year.

The second exercise plays with two trends for the priority of pumped hydropower energy storage. The first linearly changes the priority from 2020 (4.71) up to the maximum (10) at the end of the simulation (2030). The second keeps the priority of 2017 constant. Additionally, the installed capacities of VRES are expanded 2.36 times in the same period to generate more over-productive hours in the year, stimulating PHES. For all the future years, the rest of the parameters remain constant with the calibrated values of 2017.

Figure 4 reveals the impact of the growth of the variable renewables introduced in the scenario. The expansion leads to a high amount of curtailed energy, reaching around 8% of the electricity demand by 2030. The PHES-priority trend with increasing values (“increasing” label in the figure) covers higher market shares in the electricity mix. However, the consumption related to the pumping mode increases, which leads to higher generations of other dispatchable plants at certain hours to cover the increments in the load demands (CHP, dam hydropower, and fossil fuels).

Figure 4.

Evolution of some suppliers and curtailment facing a renewable growth. The brown vertical line separates the historical period (left, one known point per year) from the scenario (right, two possible points per year). “_increasing” points out variables of the scenario with increasing PHES priority (up to 10 by 2030)

The positive residual load value (PRL) in over-productive hours is zero, so the potential generation of PHES is entirely available for the negative RLD market. This is why PHES has a privileged role in the system. Consequently, the priority of the pumping mode power does not affect curtailment. The effect on curtailment comes from the installed capacity of this option.

The priorities of installed capacities are used as previously to reproduce different historical conditions (2017-2020). These conditions are assumed to justify legislative and technical dynamic processes, allowing us to estimate a range of potential benefits when promoting one or other technology over the energy transition.

This subsection analyzes the objective scenario of the Spanish NECP [47], an integrated report to meet the EU’s climate and energy targets for the period from 2021 to 2030 (EU/2018/1999). The values and correspondence across categories to adapt the HPSS model to the NECP are shown in Table 4. Priorities are not moved after 2020, running four simulations according to the historical years previously calibrated.

As can be seen in Figure 5, renewables do not go further than 61% in the share of any of the four simulations. This finding relies on the specific context, i.e., the NECP Spanish report, and it may be different from one region to another. For example, a recent study achieves a renewable penetration of 80% in the electricity mix when simulating the official Danish pathway to 2030 (objective of 70% reduction in the greenhouse gas emissions) [48]. On the 2030 Australian power system represented by De Rosa et al. [49], the indicator was estimated in the range of 21.4-36.3%.

Curtailment increases with the higher penetration of variable renewables. This is caused by the lack of sufficient flexibility in the system to integrate the excess of non-dispatchable electricity production into new demands in other energy sectors, such as heating and transport. In addition, the whole 2017 was very dry, so very productive for solar and unproductive for hydropower plants. This combination generates a less flexible hourly profile in the supply side of the power system. Most of the overproduction is therefore used by storage. The profiles of nuclear and international interconnections are exogenous, a fact that further contributes to increasing the likelihood of over-productive hours. Fortunately, although the conditions of 2017 seem to be harder than other historical years, curtailment does not so increase very much, wasting energy in 519 of the hours in 2030 (1.4 % of the VRES generation) in the worst case.

The progression of renewables in the mix reaches an upper limit. This suggests additional installations of flexibility options to further integrate the excess of electricity production and avoid carbon emissions.

Figure 5.

Percentages of renewable share in the mix and curtailment (related to VRES generation) for the NECP scenario, based on the 4 historical years

The last simulation evaluates a way to achieve the objective of a 74% renewable share in the Spanish power system. This plan assumes the following measures:

• Constant electricity demand over time.

• Reduction of the international interconnection capacity to test the Iberian case, given “the lack of a coordinated common energy policy with its neighboring systems” [50].

• Decommissioning of nuclear power plants, since they introduce inflexible generation into the model.

Finally, the optimal ratio wind/solar-PV of 1.86 is applied from the literature [51]. This ratio reduces the statistical renewable variability burden by installing infrastructure according to the regional conditions of Spain (59520 MW and 32000 MW, respectively, by 2030).

The resulting share of renewable electricity in the mix is shown in Figure 6, where 1.7 % of its generation would be curtailed in the worst conditions (2017), and which is far below 5%, a standard limit in the literature of energy planning [52],[53]. Statistically, the conditions for 2017 and 2020 would help more than those of 2018 and 2019, to reach the objective.

Figure 6.

Proposal of scenario for Spain between 2020-2030. Percentages of renewable share, percentage of curtailment (related to VRES generation), based on the assumption of 4 different historical conditions

As shown before, the generation of PHES is sensitive to changes in its priority. One way to address the uncertainty of this model parameter is through sensitivity simulation [54]. Uniform distribution from 1 to 10 is applied to the PHES priority to generate 2000 simulations. Figure 7 shows a likelihood range between 63-78% in the share of renewables in the mix, highlighting a significant influence of this technology when its priority is widely modified.

Figure 7.

Sensitivity analysis to the PHES priority from 1 to 10, with a total of 2000 simulations performed by a random uniform distribution

Hourly results reflect additional information. For instance, Figure 8 depicts the simulation with the historical conditions of 2020 at day 69 in 2030. Wind and solar contribute every hour, increasing during the daytime. In hours with low levels of renewables, thermal power plants increased the generation, whose hourly-average maximum ramp-up would increase by 1699 MW, while the ramp-down would decrease by -1319 MW.

As mentioned in the methods section, PHES runs in over-productive and deficit hours, as h1639. This situation of witnessing pumping and turbine modes of the PHES units at the same time can be found in the real operation of the power system, given the large geographical area analyzed. However, it is important to note that this result is not derived from findings with a geographical tool, but from the statistical method applied to distribute the consumption of the pumping mode (PHES_pos_cons).

Figure 8.

Hourly results from the dispatch allocation in hours 1633 to 1656 (day 69) of 2030 (proposal scenario). PHES in pumping mode is depicted in the negative axis, while PHES in turbine mode is in the positive axis