Computation of electricity marginal cost per hour with the Merit order using Rstudio
Oyeniyi Olaniyi
3/28/2020
Computation of electricity marginal cost per hour
After the second semester of my master program at TU Berlin, I thought I had understood at least the most important basics of energy engineering and their systems. I was totally not wrong. At that time, I had gotten a good grasp of major technical advantages and disadvantages of both large-scale conventional/renewable energy conversion systems and small-scale energy systems, e.g. mini-grids and swarm grids (looking forward to writing a comparison on the two soon). When I started the third semester, I began to apprehend something new, I had never thought of, a new dive into the energy ocean. I used to think all the power plants were put online always or at once to serve the energy demand specific to an individual country, for instance Germany. I did not even know anything about load duration curve or the variation of a nation’s energy demand during the year.
Using Germany as an example, I will try explaining the simplified task of the transmission system operators (TSOs) regarding power plants sequencing to meet the energy demand (in this example, per hour) and how the TSOs calculate the marginal cost of electricity in these hours. Before we can follow on the electricity marginal cost concept, the following questions will be answered, who are TSOs? What is their basic approach in a generation - load balancing? what is residual load? and, what is merit order?
Transmission System Operators (TSOs): They are entrusted with the task of transporting electricity and in some cases, natural gas using their fixed infrastructure as defined by the European commission Source. To explain a bit further, the electricity system of a country is fundalmentally divided into generation (energy suppliers), transmission (transport electricity from suppliers over high voltage lines to the distributors) and distribution (step down high voltage electricity and supply to end-users). In Germany there are four TSOs, namely, 50 Hertz, Amprion, TenneT and TransnetBW.
Generation-Demand balancing: In relation to generation-demand balancing, the TSOs with the predicted load demand for the next 24 hours (predictions are often corrected within upto 15 minutes interval) and with the weather forecast, predict the non-dispatchable renewable energy, e.g. solar and wind, for the same period. The sum of non-dispatchable renewable energy and dispatchable renewable energy, e.g. biomass is subtracted from the predicted energy demand and the remaining power demand is the residual load covered by the conventional power plants. This way priority is given to renewable energy suppliers.
Residual load = Demand load - (Non - dispatchable RES + Dispatchable RES)
library(readxl)
library(ggplot2)
library(tidyverse)## -- Attaching packages -------------------------------------------------------------------- tidyverse 1.3.0 --## v tibble 2.1.3 v dplyr 0.8.4
## v tidyr 1.0.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## v purrr 0.3.3## -- Conflicts ----------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()Total_LoadData <- suppressWarnings(read_excel("TotalLoad_DayAhead_Actual.xlsx",
sheet = "Sheet2", col_types
=c("numeric", "text", "numeric", "numeric")))
GenerationWindSolar <- suppressWarnings(read_excel("GenerationWind_and_Solar.xlsx",
sheet = "Sheet2",
col_types = c("text","numeric",
"numeric", "numeric",
"numeric")))
MinutesInterval = Total_LoadData$Time
Day = rep("12/03/2020", length(MinutesInterval))
Day[length(Day)] = "13/03/2020"
LoadDemand = Total_LoadData$`Actual Total Load`
#### Just using Solar and Wind power generation for RES consideration
RESgeneration = GenerationWindSolar$SolarWind_total
ResidualLoad = LoadDemand - RESgeneration
minutes_array = numeric(length(LoadDemand))
for (minute_int in 1:length(LoadDemand)){
if (minute_int == 1){
minutes_array[minute_int] = 15
}
else{minutes_array[minute_int] = minutes_array[minute_int-1] + 15}
}
DemandToResidualLoad = cbind.data.frame(Day, minutes_array, LoadDemand, RESgeneration,
ResidualLoad)
DemandToResidualLoad %>%
ggplot(mapping = aes(x = minutes_array,
y = LoadDemand)) +
geom_point(mapping = aes(colour = "LoadDemand")) +
geom_point(mapping = aes(y = RESgeneration, colour = "RESgeneration")) +
geom_point(mapping = aes(y = ResidualLoad, colour = "ResidualLoad")) +
ggtitle("12/03/2020", subtitle = "Data source - https://transparency.entsoe.eu") +
scale_x_continuous("Minute Interval (min)") +
scale_y_continuous("Power (MW)") +
theme(axis.title.x = element_text(face="bold", colour="black", size=14),
axis.text.x = element_text(angle=360 ,vjust=0.5, size=11, face="bold")) +
theme(axis.title.y = element_text(face="bold", colour="black", size=14),
axis.text.y = element_text(angle=360, vjust=0.5, size=11,face="bold" )) +
theme(plot.title = element_text(colour="black", size=13, face="bold", hjust = 0.5)) 
Merit order: The merit order ranks available sources (conventional power plants) of electrical generation based on their marginal cost and the cumulative amount of energy generated. This way, the TSO ensures that the plants with the lowest marginal cost are first brought online with the goal of minimizing the cost of electricity production Source. The graph below illustrates the merit order depicting the marginal cost of the power plants against the power cumulative.
Please note that not all the plants in the data source were taken for this analysis but just enough to cover residual load, especially currently functioning power plants
ESP_data = read_excel("power_plants.xlsx",
sheet = "Merit Order")
#### Data source - https://www.bundesnetzagentur.de
#### Please note that not all the plants in the data source were taken for this analysis
### but just enough to cover residual load, especially currently functioning power plants
ESP_data$Fuel = ifelse(ESP_data$Fuel == "Braunkohle",
"Brown coal", ESP_data$Fuel )
ESP_data$Fuel =ifelse(ESP_data$Fuel == "Erdgas",
"Natural gas", ESP_data$Fuel )
ESP_data$Fuel = ifelse(ESP_data$Fuel == "Steinkohle",
"Hard coal", ESP_data$Fuel )
ESP_data %>%
ggplot(mapping = aes(x = `Power cum`,
y = `Marginal cost of the plant` ,
colour = Fuel)) +
geom_col(size = 10) +
ggtitle("Merit Order", subtitle = "Data source - Bundesnetzagentur.de") +
scale_x_continuous("Power cummulative (MW)") +
scale_y_continuous("Marginal cost (€/MWh)") +
theme(axis.title.x = element_text(face="bold", colour="black", size=14),
axis.text.x = element_text(angle=360 ,vjust=0.5, size=11, face="bold")) +
theme(axis.title.y = element_text(face="bold", colour="black", size=14),
axis.text.y = element_text(angle=360, vjust=0.5, size=11,face="bold" )) +
theme(legend.text = element_text(colour="black", size = 12, face = "bold")) +
theme(plot.title = element_text(colour="black", size=13, face="bold", hjust = 0.5)) 
ESP_data %>%
rename(PlantConfiguration = Type) %>%
ggplot(mapping = aes(x = `Power cum`,
y = `Marginal cost of the plant` ,
colour = PlantConfiguration)) +
labs(title ="Merit Order with 15€/kgCO2",
subtitle = "Data source - Bundesnetzagentur.de") +
geom_col(size = 10) +
scale_x_continuous("Power cummulative (MW)") +
scale_y_continuous("Marginal cost (€/MWh)") +
theme(axis.title.x = element_text(face="bold", colour="black", size=14),
axis.text.x = element_text(angle=360 ,vjust=0.5, size=11, face="bold")) +
theme(axis.title.y = element_text(face="bold", colour="black", size=14),
axis.text.y = element_text(angle=360, vjust=0.5, size=11,face="bold" )) +
theme(legend.text = element_text(colour="black", size = 12, face = "bold")) +
theme(legend.title = element_text(colour="black", size=13, face="bold")) +
theme(plot.title = element_text(colour="black", size=13, face="bold", hjust = 0.5)) 
The two graphs are similar but the first depicts the merit order with the fuel types and the second, merit order with the different thermodynamic cycles(plant configurations).
The marginal cost of each plant is calculated considering the variable (OPEX) cost (€/MWh), fuel price (€/MWh), carbon emission (kgCO2/MWh) and carbon price (€/kgCO2) in the following formula.
Marginal cost (MC) = (OPEXvar + fuel price + (CO2emm * CO2cert))/Energy_efficiency
We can observe that the natural gas-fired power plants have higher marginal cost, majorly due to the high cost of the fuel compared to coal.
There is one more thing we observe from the graphs above, there is no RES portrayed in the merit order. In Germany, they are prioritized as explained above with the calculation of residual load. In addition, though I had included combined cycle plants in the merit order, they are usually given preference due to their supporting product (heat) needed during winter/cold season. Hence they do not necessary partake in the merit order and are often online before conventional plants of other plant configuration.
One very important point from the merit order is that the lower the marginal cost of a plant the more profitable it is since the marginal plant (last plant online to meet the required demand) sets the marginal cost of electricity. Thus running plants before the marginal plant, gain the difference of the marginal cost of electricity and their respective marginal cost of electricity production.
According to Fraunhofer ISI, the main driving factors of the merit order are the development of the fuel prices, CO2 prices and inclusion of the renewables. The first two factors we can see their influence in the marginal cost formula. The last, pushes the base load plants (coal fired power plants) towards peak load plants. This shift of coal fired power plants which are not as flexible (ramp rate) as the gas fired power plants, reduces the part load efficiencies of this coal plants, therefore increasing the marginal cost of produced electricity and increases these plants’ CO2 emissions. Despite President Donald Trump vowing to save the coal industry, the inclusion of renewables makes it difficult for them to profit on long term basis, as the market itself is making investment in coal technologies unfavorable, with coal plants becoming marginal plants more often than they ever were, hence not making much profits. In addition, inclusion of renewables affects gas plants too, as they run less and pushed up higher in the curve. To read further
Effect of four (4) times increase in CO2 certificate prices (60€/kgCO2)
ESPHighCO2 <- ESP_data <- read_excel("power_plantsHighCO2.xlsx",
sheet = "Merit Order")
ESPHighCO2$Fuel = ifelse(ESPHighCO2$Fuel == "Braunkohle",
"Brown coal", ESPHighCO2$Fuel )
ESPHighCO2$Fuel =ifelse(ESPHighCO2$Fuel == "Erdgas",
"Natural gas", ESPHighCO2$Fuel )
ESPHighCO2$Fuel = ifelse(ESPHighCO2$Fuel == "Steinkohle",
"Hard coal", ESPHighCO2$Fuel )
ESPHighCO2 %>%
ggplot(mapping = aes(x = `Power cum`,
y = `Marginal cost of the plant` ,
colour = Fuel)) +
geom_col(size = 5) +
labs(title ="Merit Order with 60€/kgCO2",
subtitle = "Data source - Bundesnetzagentur.de") +
scale_x_continuous("Power cummulative (MW)") +
scale_y_continuous("Marginal cost (€/MWh)") +
theme(axis.title.x = element_text(face="bold", colour="black", size=14),
axis.text.x = element_text(angle=360 ,vjust=0.5, size=11, face="bold")) +
theme(axis.title.y = element_text(face="bold", colour="black", size=14),
axis.text.y = element_text(angle=360, vjust=0.5, size=11,face="bold" )) +
theme(legend.text = element_text(colour="black", size = 12, face = "bold")) +
theme(plot.title = element_text(colour="black", size=13, face="bold", hjust = 0.5)) 
We can spot the shift of hard coal power plants towards the peak region and the reverse for gas power plants meaning they will run less during the year and with lower part load efficiencies. This way with the inclusion of renewables and increase in CO2 prices, some natural gas-powered plants will operate more often in the baseload region and become appealing to investors. Furthermore, a spike increase in the cost of electricity for end-users is likely to occur.
##### Convert 15 minute interval data to Hour intervals
End_hour = 24
Demand_array = numeric(length(ResidualLoad)) ####In place of appending
for (minutes in 1:length(ResidualLoad)){
for (hours in 1:End_hour){
if (hours == 1){
Demand_array = replace(Demand_array, c(1:4), mean(ResidualLoad[1:4]))
}
else{
start_minute = ((hours-1)*4) + 1
end_minute = hours*4
Demand_array= replace(Demand_array, c(start_minute:end_minute),
mean(ResidualLoad[start_minute:end_minute]))
}
}
}
ResidualLoadbyHour = numeric(End_hour)
for (hour in 1:End_hour){
ResidualLoadbyHour[hour] = Demand_array[hour*4]
}
Dayhours = 1:24
ggplot(mapping = aes(x = Dayhours, y = ResidualLoadbyHour)) +
geom_point(colour = "red") +
ggtitle("Residual load against Hours", subtitle = "Data source - https://transparency.entsoe.eu") +
scale_x_discrete("12.03.2020 00:00 - 13.03.2020 00:00") +
scale_y_continuous("Residual Load (MW)")+
theme(axis.title.x = element_text(face="bold", colour="black", size=11),
axis.text.x = element_text(angle=360 ,vjust=0.5, size=11, face="bold")) +
theme(axis.title.y = element_text(face="bold", colour="black", size=12),
axis.text.y = element_text(angle=360, vjust=0.5, size=11,face="bold" )) +
theme(plot.title = element_text(colour="black", size=11, face="bold", hjust = 0.5)) 
Electricity marginal cost
The residual load of each hour is compared to the power cummulative of the merit order, with the aim of finding the plants that are online, with respective power cummulative lower than the residual load, and offline, Power cum >=ResidualLoad.
Having the number of power plants offline, the relative part load percentages are calculated for all the plants. Then, using a polymonial model function containing an estimated relationship between the part load efficiencies and the relative part load percentages, the part load efficiencies of all the offline plants were computed with respect to the plant configuration, steam turbine (Rankine cycle), combined cycle and open gas turbine cycle.
Plants which do not fulfill the required power addition, are not considered and with reference to the plant configuration, the minimum part load percentages; for instance, gas turbines operating below 15% was set to the minimum part load percentage, 15%.
The respective marginal cost of the offline plants was divided by the corresponding part load efficiencies and the minimum marginal cost computed was the marginal cost of electricity for the residual load at the corresponding hour.
Please note the difference between the part load percentages and the part load efficiencies
names(ESP_data)## [1] "Kraftwerksnummer Bundesnetzagentur" "Plant name"
## [3] "Fuel" "Type"
## [5] "Net Power" "Marginal cost of the plant"
## [7] "Power cum"steam_part_load = c(30, 40, 50, 60,70, 80,90, 100)
steam_efficiency = c(90, 94, 97, 98,98, 98, 99, 100)
modelSteam = lm(steam_efficiency ~ poly(steam_part_load, 3))
combine_part_load = c(40, 50, 60,70, 80,90, 100)
combine_efficiency = c(83, 87, 91, 95, 97, 98, 100)
modelCombine = lm(combine_efficiency ~ poly(combine_part_load, 3))
OpenGas_part_load = c(30, 40, 50, 60,70, 80,90, 100)
OpenGas_efficiency = c(62, 72, 82, 85, 87, 93, 96,100)
modelOpenGas = lm(OpenGas_efficiency ~ poly(OpenGas_part_load, 3))
marginal_cost_array = c()
ResidualLoadCase = ResidualLoadbyHour
for (i in 1:length(ResidualLoadCase)){
plant_power = ESP_data$`Net Power`[which(ESP_data$`Power cum`>=ResidualLoadCase[i])]
Power_cumm = ESP_data$`Power cum`[which(ESP_data$`Power cum`>=ResidualLoadCase[i])]
PartLoad_array = c()
for (k in 1:length(Power_cumm)){
PartLoad = 100 - (Power_cumm[1]-ResidualLoadCase[i])*100/plant_power[k]
PartLoad_array = append(PartLoad_array, PartLoad)
}
PartLoad_array = ifelse(PartLoad_array<0, NA, PartLoad_array)
Part_load_effArray = c()
for (j in 1:length(PartLoad_array)){
plants_offline = ESP_data$Type[which(ESP_data$`Power cum`>=ResidualLoadCase[i])][j]
if (isTRUE(plants_offline == "Steam turbine")){
PartLoad_array[j] = ifelse(PartLoad_array[j]<40, 40, PartLoad_array[j])
PartLoad_dataFrameSteam = data.frame(steam_part_load = PartLoad_array[j])
Part_load_eff = predict(modelSteam, PartLoad_dataFrameSteam)[[1]]
}
if (isTRUE(plants_offline == "Open Gas turbine" )){
PartLoad_array[j] = ifelse(PartLoad_array[j]<15, 15, PartLoad_array[j])
PartLoad_dataFrameGas = data.frame(OpenGas_part_load = PartLoad_array[j])
Part_load_eff = predict(modelOpenGas, PartLoad_dataFrameGas)[[1]]
}
if (isTRUE(plants_offline == "Combined cycle")){
PartLoad_array[j] = ifelse(PartLoad_array[j]<20, 20, PartLoad_array[j])
PartLoad_dataFrameCombined = data.frame(combine_part_load = PartLoad_array[j])
Part_load_eff = predict(modelCombine, PartLoad_dataFrameCombined)[[1]]
}
Part_load_effArray = append(Part_load_effArray, Part_load_eff)
}
Part_load_eff = ifelse(Part_load_eff<0, NA, Part_load_eff)
MarginalCost = ESP_data$`Marginal cost of the plant`[which(ESP_data$`Power cum`>=ResidualLoadCase[i])]
marginal_cost = min(MarginalCost*100/Part_load_effArray, na.rm = TRUE)
marginal_cost_array = append(marginal_cost_array, marginal_cost)
}
DemandHour_df = cbind.data.frame(ResidualLoadbyHour, marginal_cost_array)
DemandHour_dfOfcourse you might notice, with higher power demand, the electricity price goes up and the correlation between both variables is 0.93. This means, as the demand tends towards its peak, plants with expensive marginal cost and good ramping rate come online, e.g. gas power plants.
This code was originally created for 10 columns, representing different test cases as presented by the Agora framework, with 8670 rows each, the hours of the year. To simplify the original task, I imported new energy load data and RES data for just 24 hours on the 12th of March, 2020 from Entsoe.
I hope you enjoyed reading this write-up, if you have suggestions or questions, please let me know. Vielen Dank