The final results results can be found from model run 1.11.2015 (token 144638929414). It is the final archived version in English. Objects were stored, so you can download the whole assessment to R in your own computer.
# This code is Op_en7237/finalresults on page [[Helsinki energy decision 2015]]
library(OpasnetUtils)
library(ggplot2)
openv.setN(0) # use medians instead of whole sampled distributions
# Download all pre-calculated inputs, e.g. building stock.
objects.latest("Op_en7237", code_name = "intermediates") # [[Helsinki energy decision 2015]]
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] oggplot
objects.latest("Op_en7392", code_name = "translate") # [[OpasnetUtils/Translate]] translate
objects.latest("Op_en7392", code_name = "dictionary") # [[OpasnetUtils/Translate]] dictionary
BSbase <- 24 # base_size for graph fonts. Use 12 if you use savefig to sav svg fils.
BS <- BSbase
saveobjects <- FALSE
savefigs <- FALSE
language <- "Finnish"
fi <- language == "Finnish"
savefig <- function(
fil,
path = "N:/YMAL/Publications/2015/Helsingin energiapäätös/Kuvat/",
sav = if(exists("savefigs")) savefigs else FALSE,
type = "svg",
height = 18,
width = 24,
units = "cm"
) {
if(sav) {
ggsave(paste(path, fil, ".", type, sep = ""), height = height, width = width, units = units)
}
}
allplants <- c(
'Biofuel heat plants',
'CHP diesel generators',
'Data center heat',
'Deep-drill heat',
'Hanasaari',
'Household air conditioning',
'Household air heat pumps',
'Household geothermal heat',
'Household solar',
'Katri Vala cooling',
'Katri Vala heat',
'Kellosaari back-up plant',
'Loviisa nuclear heat',
'Neste oil refinery heat',
'Salmisaari A&B',
'Salmisaari biofuel renovation',
'Sea heat pump',
'Sea heat pump for cooling',
'Small fuel oil heat plants',
'Small gas heat plants',
'Small-scale wood burning',
'Suvilahti power storage',
'Suvilahti solar',
'Vanhakaupunki museum',
'Vuosaari A',
'Vuosaari B',
'Vuosaari C biofuel',
'Wind mills'
)
shutdown <- allplants[!allplants %in% noshutdown]
shutdown <- c(shutdown, "None")
customdecisions <- data.frame(
Obs = NA,
Decision_maker = "Helen",
Decision = "PlantPolicy",
Option = "Custom",
Variable = c("plantParameters", "energyProcess"),
Cell = c(
paste("Plant:", paste(shutdown, collapse = ","), ";Parameter:Max,Min,Investment cost;Time:>2015", sep = ""),
""
),
Change = c("Replace", "Identity"),
Unit = NA,
Result = 0
)
decisions <- rbind(decisions, customdecisions)
if ("Hanasaari" %in% renovation) {
customdecisions <- data.frame(
Obs = NA,
Decision_maker = "Helen",
Decision = "PlantPolicy",
Option = "Custom",
Variable = rep(c("plantParameters", "energyProcess"), each = 2),
Cell = c(
"Plant:Hanasaari;Time:>=2018;Time:<=2060;Parameter:Max",
"Plant:Hanasaari;Time:>=2018;Parameter:Investment cost",
"Plant:Hanasaari;Fuel:Biofuel;Time:>=2018",
"Plant:Hanasaari;Fuel:Coal;Time:>=2018"
),
Change = rep(c("Replace", "Add"), each = 2),
Unit = NA,
Result = c(640, 100, -0.40, 0.40)
)
decisions <- rbind(decisions, customdecisions)
}
DecisionTableParser(decisions)
oprint(data.frame(
Running = c(
"These plants will be running in the custom plant policy:",
paste(allplants[!allplants %in% shutdown], collapse = ", ")
),
Shutdown = c(
"Plants that will be shut down in the custom plant policy:",
paste(shutdown, collapse = ", ")
)
))
EnergyNetworkDemand <- EvalOutput(EnergyNetworkDemand)
EnergyNetworkDemand@output <- rbind(
EnergyNetworkDemand@output,
data.frame(
Time = rep(c(2025, 2035, 2045, 2055, 2065), each = 4),
EnergySavingPolicy = rep(c("BAU", "Energy saving moderate", "Energy saving total", "WWF energy saving"), times = 5),
Temperature = "(-18,-15]",
EnergyConsumerDemandSource = "Formula",
EnergyConsumerDemandTotalSource = "Formula",
Fuel = "Cooling",
fuelSharesSource = "Formula",
EnergyNetworkDemandResult = 0,
EnergyNetworkDemandSource = "Formula"
)
)
#cat("All energy taxes are assumed zero.\n")
#objects.latest("Op_en4151", code_name = "fuelTax")
#fuelTax <- EvalOutput(fuelTax)
#result(fuelTax) <- 0
fuelUse <- EvalOutput(fuelUse)
fuelUse$Fuel <- factor(
fuelUse$Fuel, levels = c(
"Biofuel",
"Coal",
"Fuel oil",
"Gas",
"Light oil",
"Wood",
"Electricity",
"Electricity_taxed",
"Heat",
"Cooling"
), ordered = TRUE
)
DALY <- EvalOutput(DALYs)
DALYs <- unkeep(DALY, cols = c("Age", "Sex", "Population"))
DALYs <- oapply(DALYs, cols = c("Emission_site", "Emission_height", "Area"), FUN = sum)
DALYs <- DALYs[DALYs$Response == "Total mortality" , ]
EnergyNetworkCost <- EvalOutput(EnergyNetworkCost)
EnergyNetworkCost$Time <- as.numeric(as.character(EnergyNetworkCost$Time))
totalCost <- EvalOutput(totalCost)
totalCost@output$Time <- as.numeric(as.character(totalCost@output$Time))
totalCost <- unkeep(totalCost[totalCost$Time >= 2015 & totalCost$Time <=2065 , ], sources = TRUE)
# Net present value and effective annual cost
discount <- 0.03
times <- c(2015, 2065)
EAC <- EvalOutput(EAC)
if(saveobjects) {
objects.store(list = ls())
cat("All objects stored for later use:\n", paste(ls(), collapse = ", "), "\n")
}
############## POST_PROCESSING AND GRAPHS, VERSION FROM PERFERENCE ANALYSIS
cat(translate("NOTE! In all graphs and tables, the Total energy saving policy is assumed unless otherwise noted\n"))
cat(translate("Total DALYs/a by different combinations of policy options.\n"))
temp <- DALYs[as.character(DALYs$Time) %in% c("2015", "2035") & DALYs$Response == "Total mortality" , ]
oprint(
translate(oapply(temp, INDEX = c("Time", "EnergySavingPolicy", "PlantPolicy"), FUN = sum)),
caption = translate("Table 1: Total DALYs/a by different combinations of policy options."),
caption.placement = "top",
include.rownames = FALSE
)
bui <- oapply(buildings * 1E-6, cols = c("City_area", "buildingsSource"), FUN = sum)
bui <- truncateIndex(bui, cols = "Heating", bins = 4)
oggplot(bui[bui$EnergySavingPolicy == "BAU" , ], x = "Time", fill = "Heating", binwidth = 5) +
labs(
title = translate("Building stock in Helsinki by heating"),
y = translate("Floor area (M m2)")
)
savefig("Rakennuskannan koko Helsingissä")
oggplot(bui, x = "Time", fill = "Efficiency", binwidth = 5) +
{if(fi) facet_wrap(~ Energiansäästöpolitiikka) else facet_wrap(~ EnergySavingPolicy)} +
labs(
title = translate("Building stock in Helsinki by efficiency policy"),
y = translate("Floor area (M m2)")
)
oggplot(bui, x = "Time", fill = "Renovation", binwidth = 5) +
{if(fi) facet_wrap(~ Energiansäästöpolitiikka) else facet_wrap(~ EnergySavingPolicy)} +
labs(
title = translate("Building stock in Helsinki by renovation policy"),
y = translate("Floor area (M m2)")
)
oggplot(bui[bui$EnergySavingPolicy == "BAU" , ], x = "Time", fill = "Building", binwidth = 5) +
labs(
title = translate("Building stock in Helsinki"),
y = translate("Floor area (M m2)")
)
oggplot(buildings, x = "Time", fill = "Efficiency", binwidth = 5)+
{if(fi) facet_grid(Energiansäästöpolitiikka ~ Korjaukset) else facet_grid(EnergySavingPolicy ~ Renovation)} +
labs(
title = translate("Renovation of buildings by policy and efficiency"),
y = translate("Floor area (M m2)")
)
# Contains also other buildings than district heating and other energy than heating
hea <- EnergyConsumerDemandTotal * temperdays * 24 * 1E-3 # MW -> GWh
hea$Time <- as.numeric(as.character(hea$Time))
temp <- hea[hea$EnergySavingPolicy == "Energy saving total" & !hea$Consumable %in% c("District cooling", "Electric cooling") , ]
oggplot(truncateIndex(temp, cols = "Temperature", bins = 7), x = "Time", fill = "Temperature", binwidth = 5) +
{if(fi) facet_wrap(~ Hyödyke) else facet_wrap(~ Consumable)} +
labs(
title = translate("Energy consumption in all buildings"),
y = translate("Total energy (GWh /a)")
)
temp <- hea[!hea$Consumable %in% c("District cooling", "Electric cooling") , ]
oggplot(temp, x = "Time", fill = "Consumable", binwidth = 5) +
{if(fi) facet_wrap(~ Energiansäästöpolitiikka) else facet_wrap(~ EnergySavingPolicy)} +
labs(
title = translate("Energy consumption in all buildings"),
y = translate("Total energy (GWh /a)")
)
savefig("Helsingin vuotuinen energiantarve")
oggplot(hea, x = "Time", fill = "Consumable", binwidth = 5) +
{if(fi) facet_wrap(~ Energiansäästöpolitiikka) else facet_wrap(~ EnergySavingPolicy)} +
labs(
title = translate("Energy consumption in all buildings"),
y = translate("Total energy (GWh /a)")
)
hea2 <- EnergyNetworkDemand * temperdays * 24 / 1000 # MW -> GWh
hea2$Time <- as.numeric(as.character(hea2$Time))
oggplot(hea2, x = "Time", fill = "Fuel", binwidth = 5) +
{if(fi) facet_wrap(~ Energiansäästöpolitiikka) else facet_wrap(~ EnergySavingPolicy)} +
labs(
title = translate("Energy demand in the network"),
fill = translate("Consumable"),
y = translate("Total energy (GWh /a)")
)
savefig("Energiankulutus verkossa Helsingissä")
eb <-EnergyNetworkOptim[EnergyNetworkOptim$Process_variable_type == "Activity",]
eb <- eb[eb$EnergySavingPolicy == "Energy saving total" , ]
colnames(eb@output)[colnames(eb@output) == "Process_variable_name"] <- "Plant"
eb$Process_variable_type <- NULL
ebtemp <- eb[eb$Time %in% c("2035") & eb$PlantPolicy == "BAU" & eb$Temperature != "(-18,-15]" , ]
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 7)
oggplot(ebtemp, x = "Temperature", fill = "Plant", turnx = TRUE) +
labs(
title = translate("Power plant activity by temperature daily optim \nPlant policy = BAU, Year = 2035"),
x = translate("Temperature of the day"),
y = translate("Average daily activity (MW)")
)
ebtemp <- eb[eb$Time %in% c("2035") & eb$Temperature != "(-18,-15]" , ]
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 10)
oggplot(ebtemp, x = "Temperature", fill = "Plant", turnx = TRUE) +
{if(fi) facet_wrap(~ Voimalapolitiikka) else facet_wrap(~ PlantPolicy)} +
labs(
title = translate("Power plant activity by temperature daily optim in 2035"),
x = translate("Temperature of the day"),
y = translate("Average daily activity (MW)")
)
savefig("Helsingin päivittäinen kaukolämpötase")
ebtemp <- eb[eb$Time %in% c("2005") & eb$PlantPolicy == "BAU" & eb$Temperature == "(0,3]" , ]
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 10)
oggplot(ebtemp, x = "Plant", fill = "Plant", turnx = TRUE) +
{if(fi) facet_wrap(~ Lämpötila) else facet_wrap( ~ Temperature)} +
theme(axis.text.x = element_blank()) + # Turn text and adjust to right
labs(
title = translate("Power plant activity by temperature daily optim \nPlant policy = BAU, Year = 2005"),
y = translate("Average daily activity (MW)")
)
fu <- fuelUse / 3.6E+6 # From MJ/a -> GWh/a
fu <- fu[fu$EnergySavingPolicy == "Energy saving total" , ]
fu$Burner <- NULL
fu$Time <- as.numeric(as.character(fu$Time))
futemp <- fu[fu$Time %in% c("2015", "2035", "2065") & fu$PlantPolicy == "BAU" , ]
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7) * -1
oggplot(futemp, x = "Fuel", fill = "Plant", turnx = TRUE) +
{if(fi) facet_grid(Aika ~ Energiansäästöpolitiikka) else facet_grid(Time ~ EnergySavingPolicy)} +
labs(
title = translate("Energy commodity flows \n Plant policy = BAU"),
y = translate("Total annual energy (GWh/a)")
)
futemp <- fu[fu$Time %in% c("2005") & fu$PlantPolicy == "BAU" , ]
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7) * -1
oggplot(futemp, x = "Fuel", fill = "Plant", turnx = TRUE) +
labs(
title = translate("Energy commodity flows in 2005 \n Plant policy = BAU"),
y = translate("Total annual energy (GWh/a)")
)
futemp <- fu[fu$Fuel %in% c("Heat") , ]
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10) * -1
oggplot(futemp,
x = "Time", fill = "Plant", binwidth = 5) +
{if(fi) facet_wrap(~ Voimalapolitiikka) else facet_wrap(~ PlantPolicy)} +
labs(
title = translate("District heat flow"),
y = translate("Total annual energy (GWh/a)")
)
savefig("Helsingin vuotuinen kaukolämpötase")
futemp <- fu[fu$Fuel %in% c("Electricity") & fu$Plant != "Kymijoki River's plants", ]
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7) * -1
# Does not contain plants outside Helsinki: Kymijoki River's plants, a share of Olkiluoto nuclear plant.
oggplot(futemp, x = "Time", fill = "Plant", binwidth = 5) +
{if(fi) facet_wrap(~ Voimalapolitiikka) else facet_wrap(~ PlantPolicy)} +
labs(
title = translate("Electricity flow"),
y = translate("Total annual energy (GWh/a)")
)
savefig("Helsingin vuotuinen sähkötase")
emis <- truncateIndex(emissions, cols = "Emission_site", bins = 5)
emis <- emis[emis$EnergySavingPolicy == "Energy saving total" & emis$Fuel != "Electricity" , ]
levels(emis$Fuel)[levels(emis$Fuel) == "Electricity_taxed"] <- "Electricity bought"
emis$Time <- as.numeric(as.character(emis$Time))
oggplot(emis, x = "Time", fill = "Fuel", binwidth = 5) +
{if(fi) facet_grid(Saaste ~ Voimalapolitiikka, scale = "free_y") else facet_grid(Pollutant ~ PlantPolicy, scale = "free_y")} +
labs(
title = translate("Emissions from heating in Helsinki"),
y = translate("Emissions (ton /a)")
) +
scale_x_continuous(breaks = c(2000, 2050))
savefig("Helsingin energiantuotannon päästöt")
da <- DALYs[DALYs$EnergySavingPolicy == "Energy saving total" & DALYs$Fuel != "Electricity" , ]
levels(da$Fuel)[levels(da$Fuel) == "Electricity_taxed"] <- "Electricity bought"
da$Time <- as.numeric(as.character(da$Time))
oggplot(da, x = "Time", fill = "Fuel", binwidth = 5) +
{if(fi) facet_wrap(~ Voimalapolitiikka) else facet_wrap(~ PlantPolicy)} +
labs(
title = translate("Health effects of PM2.5 from heating in Helsinki"),
y = translate("Health effects (DALY /a)")
)
savefig("Helsingin energiantuotannon terveysvaikutukset")
fp <- fuelPrice[fuelPrice$Fuel %in% c(
"Biofuel",
"Coal",
"Electricity_taxed",
"Fuel oil",
"Heat",
"Light oil",
"Natural gas",
"Peat"
) , ]
fp$Time <- as.numeric(as.character(fp$Time))
levels(fp$Fuel)[levels(fp$Fuel) == "Electricity_taxed"] <- "Electricity"
ggplot(translate(fp@output), if(fi) {
aes(x = Aika, y = fuelPriceResult, colour = Polttoaine, group = Polttoaine)
} else {
aes(x = Time, y = fuelPriceResult, colour = Fuel, group = Fuel)
}) +
geom_line(size = 2)+theme_gray(base_size = BS) +
labs(
title = translate("Fuel prices (with tax)"),
y = translate("Price (Eur/MWh)")
)
savefig("Polttoaineiden verolliset hinnat")
tc <- truncateIndex(totalCost, cols = "Plant", bins = 11) / 10 * -1 # Yearly benefits (costs are negative)
tc <- tc[tc$EnergySavingPolicy == "Energy saving total" , ]
oggplot(tc, x = "Time", fill = "Cost", binwidth = 10) +
{if(fi) facet_wrap(~ Voimalapolitiikka) else facet_wrap( ~ PlantPolicy)} +
labs(
y = translate("Yearly cash flow (Meur)"),
title = translate("Total benefits and costs of energy production")
)+
scale_x_continuous(breaks = c(2000, 2020, 2040, 2060))
savefig("Energiantuotannon kokonaiskustannus Helsingissä kustannuksittain")
oggplot(tc, x = "Time", fill = "Plant", binwidth = 10) +
{if(fi) facet_wrap(~ Voimalapolitiikka) else facet_wrap(~ PlantPolicy)} +
labs(
y = translate("Yearly cash flow (Meur)"),
title = translate("Total benefits and costs of energy production")
)+
scale_x_continuous(breaks = c(2000, 2020, 2040, 2060))
savefig("Energiantuotannon kokonaiskustannus Helsingissä voimaloittain")
eac <- EAC[EAC$EnergySavingPolicy == "Energy saving total" , ] * -1
BS <- BSbase * 0.7 # Plot the next two graphs with smaller font because they are busy graphs.
eac2 <- eac[!eac$Plant %in% c(
'Household air conditioning',
'Household solar',
'Katri Vala cooling',
'Kellosaari back-up plant',
'Sea heat pump for cooling',
'Small-scale wood burning',
'Suvilahti power storage',
'Suvilahti solar',
'Vanhakaupunki museum',
'Wind mills'
) , ]
oggplot(eac2, x = "PlantPolicy", fill = "Cost", turnx = TRUE) +
{if(fi) facet_wrap(~ Voimala, scale = "free_y") else facet_wrap(~ Plant, scale = "free_y")} +
labs(
title = translate("Incomes and costs by plant"),
y = translate("Effective annual cash flow (Meur/a)")
)
savefig("Helsingin voimalaitosten kustannustehokkuus")
oggplot(eac2, x = "PlantPolicy", fill = "Cost", turnx = TRUE)+
{if(fi) facet_wrap(~ Voimala) else facet_wrap(~ Plant)} +
labs(
title = translate("Incomes and costs by plant"),
y = translate("Effective annual cash flow (Meur/a)")
)
savefig("Helsingin voimalaitosten kustannustehokkuus yhtenäisasteikolla")
BS <- BSbase
eac <- truncateIndex(eac, cols = "Plant", bins = 11)
oggplot(eac, x = "PlantPolicy", fill = "Plant", turnx = TRUE)+
labs(
title = translate("Incomes and costs by plant policy"),
y = translate("Effective annual cash flow (Meur/a)")
)
oggplot(eac, x = "PlantPolicy", fill = "Cost", turnx = TRUE)+
labs(
title = translate("Incomes and costs by plant policy"),
y = translate("Effective annual cash flow (Meur/a)")
)
savefig("Teholliset tulot ja menot energiantuotannosta Helsingissä kustannuksittain")
temp <- truncateIndex(plantParameters[plantParameters$Parameter == "Max" , ], cols = "Plant", bins = 11)
temp <- temp[temp$Time >= 2000 & temp$Time <=2070 , ]
oggplot(temp, x = "Time", fill = "Plant", binwidth = 1) +
{if(fi) facet_wrap(~ Voimalapolitiikka) else facet_wrap(~ PlantPolicy)} +
labs(
title = translate("Energy production capacity by plant policy"),
y = translate("Maximum capacity (MW)")
)
savefig("Energiantuotantokapasiteetin kehitys Helsingissä")
# odag() #Plots a directed acyclic graph of ovariables used in the model.
# This causes an internal error, so it must be the last row of the model.
Rationale
Causal diagram for the assessment.
Case-specific ovariables
Name is the name of ovariable that has case-specific rather than default content. Ident is the indentifier of the code that defines the case-specific ovariable. Token is the same as Ident but it uses a specific version of the code rather than the newest version. Latest is the code for an ovariable whose dependencies will be changed, i.e. who has the case-specific ovariable as parent. Get is the same as Latest but a specific version rather than the newest version is fetched.
## This code is Op_en7237/intermediates on page [[Helsinki energy decision 2015]]
library(OpasnetUtils)
library(ggplot2)
library(rgdal)
#library(maptools)
library(RColorBrewer)
#library(classInt)
#library(RgoogleMaps)
### Technical parameters
openv.setN(0) # use medians instead of whole sampled distributions
BS <- 24 # base_size = font size in graphs
figstofile <- FALSE
saveobjects <- TRUE
########################## Case-specific data and submodels
### Decisions
decisions <- opbase.data("Op_en7237", subset = "Decisions") # [[Helsinki energy decision 2015]]
oprint(decisions)
### Energy production in Helsinki
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] findrest
objects.latest("Op_en7311", code_name = "energyProcess") # [[Helsinki energy production]]
objects.latest("Op_en7311", code_name = "plantParameters") # [[Helsinki energy production]]
objects.latest("Op_en4151", code_name = "fuelPrice") # [[Prices of fuels in heat production]]
objects.latest("Op_en5141", code_name = "EnergyNetworkOptim") # [[Energy balance]] incl EnergyNetworkCost
### Building data in Helsinki
# Observation years must be defined for an assessment.
obstime <- Ovariable("obstime", data = data.frame(Obsyear = factor(seq(1985, 2065, 10), ordered = TRUE), Result = 1))
heatingShares <- 1 # Shares of different heating types in the building stock (already in the building data)
heating_before <- FALSE # Should heatingShares be calculated before renovate and timepoints (or after)?
efficiency_before <- TRUE # Should efficiencyShares be calculated before renovate and timepoints (or after)?
objects.latest("Op_en7115", code_name = "stockBuildings") # [[Building stock in Helsinki]]
objects.latest("Op_en7115", code_name = "changeBuildings") # [[Building stock in Helsinki]]
objects.latest("Op_en7115", code_name = "demolitionRate") # [[Building stock in Helsinki]]
objects.latest("Op_en7115", code_name = "renovationRate") # [[Building stock in Helsinki]]
objects.latest("Op_en7115", code_name = "renovationShares") # [[Building stock in Helsinki]]
objects.latest("Op_en5488", code_name = "efficiencyShares") # [[Energy use of buildings]]
objects.latest("Op_en6289", code_name = "buildingstest") # [[Building model]] # Generic building model.
## Energy use
#objects.latest("Op_en5488", code_name = "energyUseTemperature") # [[Energy use of buildings]] temperature version of energyUse
objects.latest("Op_en5488", code_name = "temperene") # [[Energy use of buildings]] temperature version of energyUse
objects.latest("Op_en5488", code_name = "nontemperene") # [[Energy use of buildings]] temperature version of energyUse
objects.latest("Op_en2959", code_name = "temperatures") # [[Outdoor air temperature in Finland]]
#objects.latest("Op_en5488", code_name = "energyUseTemperature") # [[Energy use of buildings]]
objects.latest("Op_en5488", code_name = "EnergyConsumerDemand") # [[Energy use of buildings]]
requiredName <- "Heat" # Name of the fuel type that must match for input and output
### Energy and emissions
objects.latest("Op_en7311", code_name = "emissionLocationsPerPlant") # [[Helsinki energy production]] also heatingShares
objects.latest("Op_en2791", code_name = "emissionFactors") # [[Emission factors for burning processes]]
objects.latest("Op_en2791", code_name = "emissionstest") # [[Emission factors for burning processes]]
objects.latest("Op_en7311", code_name = "fuelShares") # [[Helsinki energy production]]
objects.latest("Op_en7311", code_name = "nondynsupply") # [[Helsinki energy production]]
timelylimit <- 10000 # The largest production (MW) available at each plant (used to describe time-varying limits).
## Exposure
objects.latest("Op_en5813", code_name = "exposure") # [[Intake fractions of PM]] uses Humbert iF as default.
## Health assessment
objects.latest('Op_en2261', code_name = 'totcases') # [[Health impact assessment]] totcases and dependencies.
objects.latest('Op_en5461', code_name = 'DALYs') # [[Climate change policies and health in Kuopio]] DALYs, DW, L
population <- 623732 # Contains only the Helsinki city, i.e. assumes no exposure outside city. (Wikipedia)
# Note: the population size does NOT affect the health impact if the exposure-response function in linear.
# However, it DOES affect exposure estimates.
# DALYshortcut is for a quicker health impact calculations, as it directly uses the Helsinki conditions.
objects.latest("Op_en7237", code_name = "DALYshortcut") # [[Helsinki energy decision 2015]]
objects.latest("Op_en7379", code_name = "externalities") # [[External cost]]
objects.latest("Op_en3283", code_name = "totalCost") # [[Economic impacts]]
objects.latest("Op_en3283", code_name = "EAC") # [[Economic impacts]]
################################# Actual model and some intermediate processing.
DecisionTableParser(decisions)
# Remove previous decisions, if any.
forgetDecisions()
renovationRate <- EvalOutput(renovationRate) * 10 # Rates for 10-year periods
buildings <- EvalOutput(buildings)
buildings$EnergySavingPolicy <- factor(
buildings$EnergySavingPolicy,
levels = c("BAU", "Energy saving moderate", "Energy saving total", "WWF energy saving"),
ordered = TRUE
)
EnergyConsumerDemand <- EvalOutput(EnergyConsumerDemand)
EnergyNetworkDemand <- EvalOutput(EnergyNetworkDemand)
EnergyNetworkDemand@output <- rbind(
EnergyNetworkDemand@output,
data.frame(
Time = rep(c(2025, 2035, 2045, 2055, 2065), each = 4),
EnergySavingPolicy = rep(c("BAU", "Energy saving moderate", "Energy saving total", "WWF energy saving"), times = 5),
Temperature = "(-18,-15]",
EnergyConsumerDemandSource = "Formula",
EnergyConsumerDemandTotalSource = "Formula",
Fuel = "Cooling",
fuelSharesSource = "Formula",
EnergyNetworkDemandResult = 0,
EnergyNetworkDemandSource = "Formula"
)
)
########################### SAVE OBJECTS
if(saveobjects) {
# Clean decisions and previous results not wanted/needed by the half-model
energyProcess@output <- data.frame()
plantParameters@output <- data.frame()
EnergyNetworkOptim@output <- data.frame()
#EnergyConsumerDemand@output <- data.frame()
totcases@output <- data.frame()
DALYs@output <- data.frame()
exposure@output <- data.frame()
rm(saveobjects, dictionary) # Remove technical objects that may be updated independently of the model
objects.store(list = ls()) # Save all objects of the global environment.
cat("All objects stored for later use:\n", paste(ls(), collapse = ", "), "\n")
}
############################ Output tables and graphs
if(FALSE) {
fuelUse <- EvalOutput(fuelUse)
fuelUse$Fuel <- factor(
fuelUse$Fuel, levels = c(
"Biofuel",
"Coal",
"Fuel oil",
"Gas",
"Light oil",
"Wood",
"Electricity",
"Electricity_taxed",
"Heat",
"Cooling"
), ordered = TRUE
)
# Calculate exposure based on average iF.
exposure <- EvalOutput(exposure)
exposure <- exposure[exposure$Area == "Average" , ]
exposure <- oapply(exposure, cols = c("Plant", "Emission_site", "Emission_height", "Area"), FUN = sum)
totcases <- EvalOutput(totcases)
DALYs <- EvalOutput(DALYs)
EnergyNetworkCost <- EvalOutput(EnergyNetworkCost)
cat("Total DALYs/a by different combinations of policy options.\n")
temp <- DALYs[as.character(DALYs$Time) %in% c("2015", "2035", "2065") & DALYs$Response == "Total mortality" , ]
oprint(
oapply(temp, INDEX = c("Time", "EnergySavingPolicy", "PlantPolicy"), FUN = sum),
caption = "Table 1: Total DALYs/a by different combinations of policy options.",
caption.placement = "top",
include.rownames = FALSE
)
####################### Post-processing and graphs
bui <- oapply(buildings * 1E-6, cols = c("City_area", "buildingsSource"), FUN = sum)
bui <- truncateIndex(bui, cols = "Heating", bins = 4)
oggplot(bui[bui$EnergySavingPolicy == "BAU" , ], x = "Time", fill = "Heating", binwidth = 5) +
labs(
title = "Building stock in Helsinki by heating",
y = "Floor area (M m2)"
)
oggplot(bui, x = "Time", fill = "Efficiency", binwidth = 5) +
facet_grid(. ~ EnergySavingPolicy) +
labs(
title = "Building stock in Helsinki by efficiency policy",
y = "Floor area (M m2)"
)
oggplot(bui, x = "Time", fill = "Renovation", binwidth = 5) +
facet_grid(. ~ EnergySavingPolicy) +
labs(
title = "Building stock in Helsinki by renovation policy",
y = "Floor area (M m2)"
)
oggplot(bui[bui$EnergySavingPolicy == "BAU" , ], x = "Time", fill = "Building", binwidth = 5) +
labs(
title = "Building stock in Helsinki",
y = "Floor area (M m2)"
)
oggplot(buildings, x = "Time", fill = "Efficiency", binwidth = 5)+
facet_grid(EnergySavingPolicy ~ Renovation) +
labs(
title = "Renovation of buildings by policy and efficiency",
y = "Floor area (M m2)"
)
emissions$Time <- as.numeric(as.character(emissions$Time))
# Plot energy need and emissions
hea <- EnergyConsumerDemand * temperdays * 24 * 1E-9 # From W -> GWh /a.
hea <- hea[hea$Consumable == "Heating" , ]
hea <- oapply(hea, cols = c("City_area", "buildingsSource"), FUN = sum)
hea <- truncateIndex(hea, cols = "Heating", bins = 4)
oggplot(hea, x = "Time", fill = "Renovation", binwidth = 5) +
facet_wrap( ~ EnergySavingPolicy) +
labs(
title = "Energy used in heating in Helsinki",
x = "Time",
y = "Heating energy (GWh /a)"
)
eb <-EnergyNetworkOptim[EnergyNetworkOptim$Process_variable_type == "Activity",]
colnames(eb@output)[colnames(eb@output) == "Process_variable_name"] <- "Plant"
eb$Process_variable_type <- NULL
ebtemp <- eb[eb$Time %in% c("2035") & eb$EnergySavingPolicy == "BAU" & eb$PlantPolicy == "BAU" , ]
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 7)
oggplot(ebtemp, x = "Plant", fill = "Plant") + facet_wrap( ~ Temperature) +
theme(axis.text.x = element_blank()) + # Turn text and adjust to right
labs(
title = "Power plant activity by temperature daily optim \nPlant policy = BAU, Year = 2035",
y = "Power output daily average (MW)"
)
fu <- fuelUse
fu$Burner <- NULL
fu <- fu / 3.6E+6 # From MJ/a -> GWh/a
fu$Time <- as.numeric(as.character(fu$Time))
futemp <- fu[fu$Time %in% c("2015", "2035", "2065") & fu$PlantPolicy == "BAU" , ]
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7) * -1
oggplot(futemp, x = "Fuel", fill = "Plant", turnx = TRUE) + facet_grid(Time ~ EnergySavingPolicy) +
labs(
title = "Energy commodity flows \n Plant policy = BAU",
y = "Total annual energy (GWh/a)"
)
futemp <- fu[fu$Fuel %in% c("Heat") , ]
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10) * -1
oggplot(futemp, x = "Time", fill = "Plant", binwidth = 5) + facet_grid(EnergySavingPolicy ~ PlantPolicy) +
labs(
title = "District heat flow",
y = "Total annual energy (GWh/a)"
)
futemp <- fu[fu$Fuel %in% c("Electricity") , ]
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10) * -1
oggplot(futemp, x = "Time", fill = "Plant", binwidth = 5) + facet_grid(EnergySavingPolicy ~ PlantPolicy) +
labs(
title = "Electricity flow",
y = "Total annual energy (GWh/a)"
)
oggplot(emissions[emissions$EnergySavingPolicy == "BAU" , ], x = "Time", fill = "Fuel", binwidth = 5) +
facet_grid(Pollutant ~ PlantPolicy, scale = "free_y") +
labs(
title = "Emissions from heating in Helsinki",
x = "Time (Energy saving policy = BAU)",
y = "Emissions (ton /a)"
)
oggplot(exposure[exposure$Exposure_agent == "PM2.5" , ], x = "Time", fill = "Fuel", binwidth = 5) +
facet_grid(EnergySavingPolicy ~ PlantPolicy) +
labs(
title = "Exposure to PM2.5 from heating in Helsinki",
x = "Time",
y = "Average PM2.5 (µg/m3)"
)
oggplot(DALYs, x = "Time", fill = "Fuel", binwidth = 5) + facet_grid(EnergySavingPolicy ~ PlantPolicy) +
labs(
title = "Health effects in DALYs of PM2.5 from heating in Helsinki",
y = "Health effects (DALY /a)"
)
} # END if(FALSE)
Building stock in Helsinki
Question
What is the building stock in Helsinki and its projected future?
Answer
Current building stock in Helsinki by heating type.
Projected building stock based on 2015 data and urban plans.
## This code is Op_en7115/ on page [[Building stock in Helsinki]].
library(OpasnetUtils)
library(ggplot2)
objects.latest("Op_en7115", code_name = "stockBuildings")
stockBuildings <- EvalOutput(stockBuildings)
ggplot(stockBuildings@output, aes(x = Time, weight = stockBuildingsResult, fill = Heating)) +
geom_bar(binwidth = 5) + theme_gray(base_size = 24) +
labs(
title = "Current building stock (floor area) by heating type \n and year of construction",
x = "Construction year",
y = expression("Floor area ( "*m^2*" )"))
ggplot(stockBuildings@output, aes(x = Time, weight = stockBuildingsResult, fill = Building)) +
geom_bar(binwidth = 5) + theme_gray(base_size = 24) +
labs(
title = "Current building stock (floor area) by building type \n and year of construction",
x = "Construction year",
y = expression("Floor area ( "*m^2*" )"))
Rationale
This part contains the data needed for calculations about the building stock in Helsinki. It shows the different building and heating types in Helsinki, and how much and what kind of renovations are done for the existing building stock in a year, including how much and how old building stock is demolished. This data is used in further calculations in the model.
There is also some other important data that wasn't used in the model's calculations. These include more accurate renovation statistics for residential buildings, U-value changes for renovations and thermal transmittance of different parts of residential buildings. This data is found under Data not used.
# This is code Op_en/7115 on page [[Building stock in Helsinki]]
library(OpasnetUtils)
library(plotly)
objects.latest("Op_en6007", code_name = "miscellaneous") # [[OpasnetUtils/Drafts]] truncateIndex
objects.latest("Op_en6007", code_name = "hnh2035") # [[OpasnetUtils/Drafts]] pushIndicatorGraph
objects.latest("Op_en7237", code_name = "intermediates") # [[Helsinki energy decision 2015]] buildings etc.
buildings <- EvalOutput(buildings)
tmp <- truncateIndex(buildings,"Building",6)
colnames(tmp@output)[colnames(tmp@output)=="EnergySavingPolicy"] <- "Scenario"
#> unique(tmp$Scenario)
#[1] BAU Energy saving moderate Energy saving total
#[4] WWF energy saving
levels(tmp$Scenario) <- c("BAU","NA","tavoite","NA")
tmp <- tmp[tmp$Scenario!="NA",]
p_buildings_b <- plot_ly(
oapply(tmp[tmp$Scenario=="BAU",], c("Building","Time"),sum)@output,
x = ~Time, y = ~buildingsResult, color = ~Building,
type = 'scatter', mode = 'lines') %>%
layout(
title="Rakennusala talotyypeittäin",
xaxis=list(title="Vuosi"),
yaxis=list(title="Rakennusala (m2)")
)
p_buildings_h <- plot_ly(
oapply(tmp[tmp$Scenario=="BAU",], c("Heating","Time"),sum)@output,
x = ~Time, y = ~buildingsResult, color = ~Heating,
type = 'scatter', mode = 'lines') %>%
layout(
title="Rakennusala lämmitysmuodoittain",
xaxis=list(title="Vuosi"),
yaxis=list(title="Rakennusala (m2)")
)
pushIndicatorGraph(p_buildings_b, "https://hnh.teamy.fi/v1/indicator/12/")
pushIndicatorGraph(p_buildings_h, "https://hnh.teamy.fi/v1/indicator/13/")
Building stock
These tables are based on FACTA database classifications and their interpretation for assessments.
This data is used for modelling. The data is large and can be seen from the Opasnet Base. Technical parts on this page are hidden for readability. Building types should match Energy use of buildings#Baseline energy consumption.
Terveydenhoidon erityislaitokset (mm. kuntoutuslaitokset)
Buildings for institutional care
89
Vapaa-ajan asunnot
Free-time residential buildings
90
Viljankuivaamot ja viljan säilytysrakennukset, siilot
Warehouses
----#: . Viimeiset 12 riviä (ilman numeroa) ovat tyyppejä jotka ovat datassa (tai ainakin vanhassa taulukossa) mutta puuttuvat Sonjan luokittelusta. --Jouni (talk) 12:57, 24 August 2015 (UTC) (type: truth; paradigms: science: comment)
For residential buildings (classes A and B) the classification is kept more detailed than for other buildings. This is because residential buildings are the biggest energy consumers in Helsinki and different classes of residential buildings are examined separately.
Reference for the classification: http://www.stat.fi/meta/luokitukset/rakennus/001-1994/koko_luokitus.html
library(OpasnetUtils)
library(ggplot2)
# [[Building stock in Helsinki]], building stock, locations by city area (in A Finnish coordinate system)
#stockBuildings <- Ovariable("stockBuildings", ddata = "Op_en7115.stock_details")
#colnames(stockBuildings@data)[colnames(stockBuildings@data) == "Built"] <- "Time"
#colnames(stockBuildings@data)[colnames(stockBuildings@data) == "Postal code"] <- "City_area"
# [[Building stock in Helsinki]]
dat <- opbase.data("Op_en7115.stock_details")[ , c(
# "Rakennus ID",
"Sijainti",
"Valmistumisaika",
# "Julkisivumateriaali",
"Käyttötarkoitus",
# "Lämmitystapa",
"Polttoaine",
# "Rakennusaine",
# "Varusteena koneellinen ilmanvaihto",
# "Perusparannus",
# "Kunta rakennuttajana",
# "Energiatehokkuusluokka",
# "Varusteena aurinkopaneeli",
"Tilavuus",
"Kokonaisala",
"Result" # Kerrosala m2
)]
colnames(dat) <- c("City_area", "Time", "Building types in Facta", "Heating types in Facta", "Tilavuus", "Kokonaisala", "Kerrosala")
dat$Time <- as.numeric(substring(dat$Time, nchar(as.character(dat$Time)) - 3))
#dat <- dat[dat$Time != 2015 , ] # This is used to compare numbers to 2014 statistics.
dat$Time <- as.numeric(as.character((cut(dat$Time, breaks = c(0, 1885 + 0:26*5), labels = as.character(1885 + 0:26*5)))))
dat$Tilavuus <- as.numeric(as.character(dat$Tilavuus))
dat$Kokonaisala <- as.numeric(as.character(dat$Kokonaisala))
dat$Kerrosala <- as.numeric(as.character(dat$Kerrosala))
build <- tidy(opbase.data("Op_en7115.building_types"))
colnames(build)[colnames(build) == "Result"] <- "Building"
heat <- tidy(opbase.data("Op_en7115.heating_types"))
colnames(heat)[colnames(heat) == "Result"] <- "Heating"
######################
# Korjaus
########################
temp <- as.character(heat$Heating)
temp[temp == "District heating"] <- "District"
temp[temp == "Light oil"] <- "Oil"
temp[temp == "Fuel oil"] <- "Oil"
heat$Heating <- temp
########################################
dat <- merge(merge(dat, build), heat)#[c("City_area", "Time", "Building", "Heating", "stockBuildingsResult")]
dat$Kerrosala[is.na(dat$Kerrosala)] <- dat$Kokonaisala[is.na(dat$Kerrosala)] * 0.8 # If floor area is missing, estimate from total area.
cat("Kerrosala ilman 2015 (m^2)\n")
oprint(aggregate(dat["Kerrosala"], by = dat["Building"], FUN = sum, na.rm = TRUE))
cat("Kokonaisala ilman 2015 (m^2)\n")
oprint(aggregate(dat["Kokonaisala"], by = dat["Building"], FUN = sum, na.rm = TRUE))
cat("Tilavuus ilman 2015 (m^3)\n")
oprint(aggregate(dat["Tilavuus"], by = dat["Building"], FUN = sum, na.rm = TRUE))
temp <- aggregate(dat["Kerrosala"], by = dat[c("Time", "Building", "Heating")], FUN =sum, na.rm = TRUE)
colnames(temp)[colnames(temp) == "Kerrosala"] <- "stockBuildingsResult"
stockBuildings <- Ovariable("stockBuildings", data = temp)
objects.store(stockBuildings)
cat("Ovariable stockBuildings stored.\n")
Construction and demolition
It is assumed that construction occurs at a constant rate so that there is an increase of 42% in 2050 compared to 2013. Energy efficiency comes from Energy use of buildings.
# This code is Op_en7115/demolitionRate on page [[Building stock in Helsinki]]
library(OpasnetUtils)
demolitionRate <- Ovariable('demolitionRate',
dependencies = data.frame(Name = "dummy"),
formula = function(...) {
temp <- tidy(opbase.data('Op_en7115', subset = 'Demolition rate'))
temp$Age <- round(as.numeric(as.character(temp$Age)))
out <- as.data.frame(approx(
temp$Age,
temp$Result,
n = (max(temp$Age) - min(temp$Age) + 1),
method = "constant"
))
colnames(out) <- c("Age", "demolitionRateResult")
out$demolitionRateResult <- out$demolitionRateResult / 100 * 10 # For ten-year intervals
out <- Ovariable("demolitionRate", output = out, marginal = c(TRUE, FALSE))
return(out)
}
)
objects.store(demolitionRate)
cat("Object demolitionRate stored.\n")
Heating type conversion
The fraction of heating types in the building stock reflects the situation at the moment of construction and not currently. The heating type conversion corrects this by changing a fraction of heating methods to a different one at different timepoints. Cumulative fraction, other timepoints will be interpolated.
library(OpasnetUtils)
renovationShares <- Ovariable("renovationShares",
dependencies = data.frame(Name = "dummy"),
formula = function(...) {
out <- Ovariable("raw", ddata = 'Op_en7115', subset = 'Popularity of renovation types')
out <- findrest((out), cols = "Renovation", total = 100) / 100
renovationyear <- Ovariable("renovationyear", data = data.frame(
Obsyear = factor(c(2015, 2025, 2035, 2045, 2055, 2065)),
Result = 1
))
out <- out * renovationyear # renovation shares repeated for every potential renovation year.
out@output$Renovation <- factor(out@output$Renovation, levels = c(
"None",
"General",
"Windows",
"Technical systems",
"Sheath reform"
), ordered = TRUE)
return(out)
}
)
objects.store(
renovationShares # Fraction of renovation type when renovation is done.
)
cat("Objects renovationShares stored.\n")
Building model
Question
How to estimate the size of the building stock of a city, including heating properties, renovations etc? The situation is followed over time, and different policies can be implemented.
Answer
Causal diagram of the building model. The actual model is up to the yellow node Building stock, and the rest is an example how the result can be used in models downstream.
The building model follows the development of a city's or area's building stock over time. The output of the model is the floor area (or volume, depending on the input data) of the building stock of a city at specified timepoints, classified by energy efficiency, heating type, and optionally by other case-specific characteristics. The model functions as part of Opasnet's modeling environment and it is coded using R. It uses specific R objects called ovariables. The model can also be downloaded and run on one's own computer.
The model is given data about the building stock of a certain city or area during a certain period of time. The data can be described with very different levels of precision depending on the situation and what kind of information is needed. Some kind of data on the energy efficiency and heating type is necessary, but even rough estimates suffice. Then again, if there is sufficient data, the model can analyse even individual buildings.
In addition to that, the model can describe changes in the building stock, i.e. construction of new buildings and demolishing
of old ones. Data on the heating- and energy efficiencies of new and demolished buildings is required at the same level of precision as that of other buildings. This data is used to calculate how construction and demolishing change the building stock's size and heating types.
The model takes into account the energy renovation of existing buildings. They are analysed using two variables:
firstly, what fraction of the building stock is energy renovated yearly and secondly, what type of renovation it is.
This information, too, can be rough or precise and detailed. It can describe the whole building stock with a single number or be
specific data on the time, the building's age, use or other background information.
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle B_{t,h,e,r} = \int\int (Bs_{c,t,a} Hs_h Es_e + Bc_{c,h,e,t,a}) Rr_a Rs_{r,t} O)\mathrm{d}c \mathrm{d}a}
B = buildings, floor area of buildings in specified groups
Bs = stockBuildings, floor area of the current buildings
Bc = changeBuildings, floor area of constructed and demolished (as negative areas) buildings
Hs = heatingShares, fractions of different heating types in a group of buildings
Es = efficiencyShares, fractions of different efficiency classes in a group of buildings
Rr = renovationRate, fraction of buildings renovated per year
Rs = renovationShares, fractions of different renovation types performed when buildings are renovated
O = obstime, timepoints for which the building stock is calculated.
Indices required (also other indices are possible)
t = Obsyear, time of observation. This is renamed Time on the output data.
c = Construction year (the index is named 'Time' in the input data), time when the building was built.
a = Age, age of building at a timepoint. This is calculated as a = t - b.
h = Heating, primary heating type of a building
e = Efficiency, efficiency class of building when built
r = Renovation, type of renovation done to a non-renovated building (currently, you can only renovate a building once)
The model is iterative across the Obsyear index so that renovations performed at one timepoint are inherited to the next timepoint, and that situation is the starting point for renovations in that timepoint.
Amount of building stock (typically in floor-m2) at given timepoints.
Required indices: Time (time the building was built. If not known, present year can be used for all buildings.) Typical indices: City_area, Building (building type)
You must give either stockBuildings, heatingShares, and efficiencyShares or changeBuildings or both. For missing data, use 0.
heatingShares (case-specific data from the user)
Fractions of heating types. Should sum up to 1 within each group defined by optional indices.
Required indices: Heating. Typical indices: Time, Building
If no data, use 1 as a placeholder.
efficiencyShares (case-specific data from the user)
Fraction of energy efficiency types. Should sum up to 1 for each group defined by other indices.
Required indices: Efficiency. Typical indices: Time, Building.
If no data, use 1 as default.
changeBuildings (case-specific data from the user)
Construction or demolition rate as floor-m2 at given timepoints.
Required indices: Obsyear, Time, Efficiency, Heating. If both stockBuildings and changeBuildings are used, changeBuildings should have all indices in stockBuildings, heatingShares, and efficiencyShares. Typical indices: Building, City_area.
If the data is only in stockBuildings, use 0 here.
renovationShares (case-specific data from the user)
Fraction of renovation types when renovation is done. Should sum to 1 for each group defined by other indices.
Required indices: Renovation, Obsyear. Obsyear is the time when the renovation is done
### This code is Op_en6289/buildingstest on page [[Building model]].
library(OpasnetUtils)
########## Calculate all building events (constructions, demolitions, renovations)
buildings <- Ovariable("buildings",
dependencies = data.frame(Name = c(
"stockBuildings",
"changeBuildings",
"heatingShares", # This can be indexed by building year (Time) or Observation year (Obsyear)
"efficiencyShares", # Same here
"renovationShares",
"renovationRate", # Fraction of buildings renovated between timepoints
"demolitionRate", # Fraction of buildings demolished between timepoints
#"heatTypeConversion", # Is dependent on buil, which defined in the formula (must handle withihn formula). Fraction of heating types converted to other
"obstime" # In the new version obstime is ovariable with column Obsyear.
)),
formula = function(...) {
stock <- stockBuildings * obstime
if( ! "Efficiency" %in% colnames(stock@output)) { # Add those shares that are missing
stock <- stock * efficiencyShares
}
if( ! "Heating" %in% colnames(stock@output)) {
stock <- stock * heatingShares
}
stock2 <- changeBuildings * obstime
if( ! "Efficiency" %in% colnames(stock2@output)) { # Add those shares that are missing
stock2 <- stock2 * efficiencyShares
}
if( ! "Heating" %in% colnames(stock2@output)) {
stock2 <- stock2 * heatingShares
}
stock2@output$Time <- as.factor(stock2@output$Time) # Also stock@output$Time is factor
buil <- OpasnetUtils::combine(stock, stock2) # Buildings in the whole timeline
buil <- unkeep(buil, sources = TRUE)
buil@output$Age <- as.numeric(as.character(buil@output$Obsyear)) - as.numeric(as.character(buil@output$Time))
buil@marginal <- c(buil@marginal, TRUE)
buil@output <- buil@output[buil@output$Age >= 0 , ]
# Note that if stockBuildings and changeBuildings have different marginals (typically policies),
# there will be NA in these indices. This will be corrected below.
for(i in colnames(buil@output)[buil@marginal]) {
if(any(is.na(buil@output[[i]]))) {
buil@output <- fillna(buil@output, i)
print(paste("Column", i, "treated with fillna (difference between stockBuildings and changeBuildings)."))
}
}
# heatTypeConversion dependent on buil, which defined above, and hence must be handled here.
# Scope is messed up due to the dependency and the following produces error when run within EvalOutput of buildings
#objects.latest("Op_en7115", code_name = "heatTypeConversion")
#heatTypeConversion <- EvalOutput(heatTypeConversion)
heatTypeConversion <- function(){
dat <- opbase.data("Op_en7115", subset = "Yearly_heating_converted_factor")
colnames(dat)[colnames(dat) == "Time"] <- "Obsyear"
dat$Obs <- NULL
out <- data.frame()
temp <- unique(dat[c("Heating_from", "Heating_to")])
for (i in 1:nrow(temp)) {
onetype <- merge(temp[i,], dat)
tempout <- merge(obstime@output, onetype, all.x = TRUE)[c("Obsyear","Result")]
tempout <- merge(tempout, temp[i,])
for (j in (1:nrow(tempout))[is.na(tempout$Result)]) {
a <- onetype$Obsyear[which.min(abs(as.numeric(as.character(onetype$Obsyear)) - as.numeric(as.character(obstime$Obsyear[j]))))]
tempout$Result[j] <- onetype$Result[a]
}
out <- rbind(out, tempout)
}
out <- Ovariable(output = out, marginal = colnames(out) != "Result")
colnames(out@output)[colnames(out@output) == "Heating_from"] <- "Heating"
out <- buil * out
out1 <- out
out1$Result <- - out1$Result
out1$Heating_to <- NULL
out$Heating <- out$Heating_to
out$Heating_to <- NULL
out@output <- rbind(out1@output, out@output)
heatTypeConversion <- out
return(out)
}
buil <- OpasnetUtils::combine(buil, heatTypeConversion())
temp1 <- merge(renovationRate, unique(buil@output["Age"])) # Avoid redundant calculations
temp1@name <- renovationRate@name
temp2 <- renovationShares
temp2@output <- temp2@output[temp2@output$Renovation == "None" , ]
renovate <- (1 - temp1) * (1 - temp2) # Assumes that data has row for Renovation: None = 0
renovate <- OpasnetUtils::combine(renovate, temp1 * renovationShares)
out <- data.frame()
prevreno <- unkeep(buil * renovate * 0, sources = TRUE) # Previously renovated buildings for the first time point
prevreno@output <- prevreno@output[prevreno@output$Renovation != "None" , ][1,] # Just take one renovated example row
marginals <- colnames(prevreno@output)[prevreno@marginal]
demolitionRate <- unkeep(demolitionRate, sources = TRUE) # Would otherwise cause trouble in rbind.
for(i in obstime@output$Obsyear) { # Accumulate the building stock
# Take the building stock
temp <- buil
temp@output <- temp@output[temp@output$Obsyear == i , ] # Take buildings of year i
prevreno@output$Obsyear <- i # update the observation year of the previous renovation
# Remove from the stock buildings that are demolished
temp <- temp * (1 - demolitionRate)
# Subtrack from the building stock buildings that are already renovated.
temp <- OpasnetUtils::combine(temp, -1 * oapply(prevreno, cols = "Renovation", FUN = sum)) # Subtract previously renovated
reno <- unkeep(temp * renovate, sources = TRUE) # Renovate the current non-renovated stock
# If renovate does not match with this timepoint, mark all buildings here as unrenovated.
# Is there new renovation in this time point?
if(sum(result(reno)) == 0) {
reno <- temp * Ovariable(
output = data.frame(Renovation = "None", Result = 1),
marginal = c(TRUE, FALSE)
)
}
newreno <- OpasnetUtils::combine(prevreno, reno) # Renovated buildings in this time point.
# Take the previous timepoints, and previously renovated and now possibly renovated together.
out <- rbind(newreno@output, out) # Add previously and now renovated together
prevreno <- unkeep(newreno, sources = TRUE)
prevreno@output <- prevreno@output[prevreno@output$Renovation != "None" , ]
}
out <- Ovariable(output = out, marginal = colnames(out) %in% marginals)
# Note that if buil and renovate have different marginals (typically policies),
# there will be NA in these indices. This will be corrected below.
for(i in colnames(out@output)[out@marginal]) {
if(any(is.na(out@output[[i]]))) {
out@output <- fillna(out@output, i)
print(paste("Column", i, "treated with fillna (difference between buil and renovate)."))
}
}
out <- oapply(unkeep(out, sources = TRUE), cols = c("Age", "Time"), FUN = sum)
colnames(out@output)[colnames(out@output) == "Obsyear"] <- "Time" # From now on, Time means the time of observation
out@output$Time <- as.numeric(as.character(out@output$Time))
return(out)
}
)
objects.store(buildings)
cat("Saved ovariable buildings\n")
Energy use of buildings
Question
How to model the use of energy of buildings based on either annual consumption per floor area, or energy efficiency per floor area per indoor-outdoor temperature difference?
Answer
Example of consumer energy demand calculations: energy need in Helsinki from the assessment Helsinki energy decision 2015.
An example code for fetching and using the variables.
## This code is Op_en5488/ on page [[Energy use of buildings]]
library(OpasnetUtils)
objects.latest("Op_en5488", code_name = "EnergyConsumerDemand") # [[Energy use in buildings]]
cat("These ovariables must be obtained from somewhere to use this method:\n")
oprint(EnergyConsumerDemand@dependencies)
Rationale
Input
Variables needed to calculate the EnergyConsumerDemand. Note that there are several different methods available, and temperature data is not needed in an annual energy version.
The code below assumes energy consumption factors relative to floor area (W /m2 /K). Local temperature data must be given in either individual or aggregated way. Individual way has temperature data for all timepoints (e.g. days or hours) of the given year, and heatingTime = 1. Aggregated way has a specific Temperature index (e.g. very cold, cold, cool etc) in both ovariables temperature and heatingTime. The ovariable temperature tells what is the actual temperature when it is "very cold", and heatingTime tells how many hours it is "very cold" during the year.
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Q_{e,r,t} = \sum_b (B_{b,e,r} U_b (17 - T_t) E_e R_r + W_b) t_t,}
where
Q = Energy used for heating and cooling (kWh /a)
B = floor area of a building stock indexed by renovation and efficiency (m2)
U = energy consumption factor per floor area for a building type (W /m2 /K)
T = temperature outside (assumes that no heating is needed if outside temperature is 17 degrees Celsius)
E = relative efficiency of a building stock based on energy class when built (no unit)
R = relative efficiency of a building stock based on energy class after renovated (no unit)
W = heating need of hot water (W)
t = time spent in a particular outdoor temperature (h /a)
### This code is Op_en5488/EnergyConsumerDemand on page [[Energy use of buildings]].
library(OpasnetUtils)
### EnergyUse of a given building stock when U values are available.
EnergyConsumerDemand <- Ovariable("EnergyConsumerDemand",
dependencies = data.frame(
Name = c(
"temperene",
"efficiencyRatio",
"renovationRatio",
"nontemperene",
"buildings",
"temperatures",
"temperdays" # fetched here but it is only calculated after energyBalance optimisation
),
Ident = c(
"Op_en5488/temperene",
"Op_en5488/efficiencyRatio",
"Op_en5488/renovationRatio",
"Op_en5488/nontemperene",
NA,
NA,
NA
)
),
formula = function(...) {
#NOTE! IF OTHER INPUTS ARE INDEXED BY Building, THIS MUST BE CHANGED!
#NOTE2! buildings must have index Heating.
#NOTE3! Heating, cooling and Other should probably be in different ovariables instead of being here:
# multiple variables necessitates explicit handling in other dependent variables, which is bad.
buil <- oapply(buildings, cols = "Building", FUN = sum)
# Remove Fuel column as unnecessary and disruptive
temperene <- temperene[,!colnames(temperene@output) %in% "Fuel"]
nontemperene <- nontemperene[,!colnames(nontemperene@output) %in% "Fuel"]
# We assume that energy efficiency of building structures also reduces cooling need.
tempe <- temperene * efficiencyRatio * renovationRatio
# We assume that 17 °C is thermoneutral with no heating.
heat <- buil * (17 - temperatures) * tempe[tempe$Consumable == "Heating", ]#!colnames(tempe@output) %in% "Fuel"]
# Only positive heating is considered (i.e., cooling is not here)
result(heat) <- pmax(0, result(heat))
#heat@output <- heat@output[heat@output$Consumable == "Heating" , ]
# We assume that 24 °C is thermoneutral with no cooling.
cool <- buil * (temperatures - 24) * tempe[tempe$Consumable %in% c("Cooling", "District cooling", "Electric cooling"),]
# Only positive cooling is considered (i.e., heating is not here)
result(cool) <- pmax(0, result(cool))
#cool@output <- cool@output[cool@output$Consumable %in% c("Cooling", "District cooling", "Electric cooling") , ]
#cool@output$Heating <- "Not heating"
other <- buil * nontemperene
#other@output$Heating <- "Not heating"
out <- OpasnetUtils::combine(heat, cool, other)
out <- unkeep(out, prevresults = TRUE, sources = TRUE)
out@output <- fillna(out@output, marginals = "Temperature")
#out[1:10]@output
#nrow(out@output)
return(out)
}
)
objects.store(EnergyConsumerDemand)
cat("Ovariable EnergyConsumerDemand stored.\n")
Baseline energy consumption
Heat reflects the energy need for heating in situations where the outdoor temperature is below 17 °C. Cooling reflects the cooling need (measured as thermal energy, not electricity!) in situations where the outdoor temperature is above 24 °C. This is not a U value, because it is about energy use per floor area, not about heat loss through building structures per m2. For estimating temperene, we take the total energy consumption in Helsinki and divide that with the total floor area and average temperature difference, see Helsinki energy consumption#U values based on overall data.
## This is code Op_en5488/temperene on page [[Energy use of buildings]]
library(OpasnetUtils)
temperene <- Ovariable("temperene", ddata = "Op_en5488", subset = "Energy use per area and temperature difference")
objects.store(temperene)
cat("Ovariable temperene stored.\n")
Temperature-independent energy consumption per floor area.
## This is code Op_en5488/nontemperene on page [[Energy use of buildings]]
library(OpasnetUtils)
nontemperene <- Ovariable("nontemperene", ddata = "Op_en5488", subset = "Temperature-independent energy use per area")
objects.store(nontemperene)
cat("Ovariable nontemperene stored.\n")
Energy efficiency in heating
What is the relative energy consumption of different efficiency classes compared with Old? This table tells that with some background information about heat (in kWh/m2/a), electricity, and water consumption.
# This code is Op_en5488/efficiencyRatio on page [[Energy use of buildings]]
library(OpasnetUtils)
efficiencyRatio <- Ovariable(
name = 'efficiencyRatio',
ddata = 'Op_en5488',
subset = 'Energy use by energy class of building'
)
objects.store(efficiencyRatio)
cat("Object efficiencyRatio initiated!\n")
# This code is Op_en5488/renovationRatio on page [[Energy use of buildings]]
library(OpasnetUtils)
renovationRatio <- Ovariable(
name = 'renovationRatio',
ddata = 'Op_en5488',
subset = 'Energy saving potential of different renovations'
)
renovationRatio@data$Renovation <- factor(
renovationRatio@data$Renovation,
levels = c("None", "General", "Windows", "Technical systems", "Sheath reform"),
ordered = TRUE
)
objects.store(renovationRatio)
cat("Object renovationRatio initiated!\n")
Helsinki energy consumption
Question
How much energy is consumed and to what purposes in Helsinki?
Answer
Energy consumed for heating, cooling, hot water, and consumer electricity in Helsinki. Note: the future consumption is based on Energy saving total scenario from the assessment Helsinki energy decision 2015, not business-as-usual.
The total heat consumption by district-heated buildings is 6921.65 GWh in 2013 (see below). We can derive the total energy efficiency value expressed as W /m2 /K for floor area and temperature difference between indoors and outdoors. The typical energy efficiency calculations (using the so called U value) assume that outdoor 17 °C is thermoneutral and lower values require heating. The total floor area of district-heated buildings is 38990000 m2 in 2015 according to the Helsinki energy decision 2015model. The annual average temperature in Helsinki is 4.8 °C [4] and during heating season Sep-May 1.4 C (Opasnet data). Therefore the energy efficiency value (approximate U value) is
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle U = \frac{6921.65 GWh/a /(24 h / d \times 365 d/a)}{38990000 m^2 (17 K - 4.8 K)} = 1.661 \frac{W}{m^2 K}}
Energy consumption statistics
Total energy consumption in Helsinki in 2013 (GWh) [1]
Adjusted for temperature
Not adjusted for temperature
District heating
6921.65
6461.00
Separate heating
303.89
284.01
Electric heating
339.23
316.65
Consumer electricity
3988.10
3988.10
Private cars
1294.06
1294.06
Other road traffic
794.33
794.33
Trains
111.16
111.16
Ships
432.12
432.12
Industry and machinery
147.60
147.60
Total
14332.14
13829.03
Helsinki energy production
Question
What is the amount of energy produced (including distributed production) in Helsinki? Where is it produced (-> emissions)? Which processes are used in its production?
## This is code Op_en7311/answer on page [[Helsinki energy production]]
library(OpasnetUtils)
objects.latest("Op_en7311", code_name = "energyProcess") # [[Helsinki energy production]]
objects.latest("Op_en7311", code_name = "plantParameters") # [[Helsinki energy production]]
objects.latest("Op_en7311", code_name = "nondynsupply") # [[Helsinki energy production]]
objects.latest("Op_en7311", code_name = "fuelShares") # [[Helsinki energy production]]
objects.latest("Op_en7311", code_name = "emissionLocationsPerPlant") # [[Helsinki energy production]]
#oprint(head(EvalOutput(energyProcess)@output))
#oprint(head(EvalOutput(plantParameters)@output))
#oprint(head(EvalOutput(nondynsupply)@output))
#oprint(head(EvalOutput(fuelShares)@output))
#oprint(head(EvalOutput(emissionLocationsPerPlant)@output))
cat("ovariables energyProcess, plantParameters, nondynsupply, fuelShares, and emissionLocationsPerPlant successfully fetched.\n")
Rationale
This page contains data about the heat plants in Helsinki. It tells, how much and what type of energy a plant produces per unit of fuel, how much the plants cost and the locations of the power plant emissions. This data is then further used in the model.
The efficiency of heat pumps is largely dependent on outside air temperature, it's feasible for a household air heat pump to reach COP 5 at 10 °C and COP 1.5 at -25 °C.
7
Household air conditioning
None
0
-0.7 - -0.2
0
1
0
0
0
0
8
Household geothermal heat
None
0
-0.36 - -0.31
1
0
0
0
0
0
Motiva 2014
9
Katri Vala cooling
None
0
-0.36 - -0.31
0
1
0
0
0
0
District cooling produced by absorption (?) heat pumps. Same as heat pumps for heating, Motiva 2014.
10
Katri Vala heat
None
0
-0.36 - -0.31
1
0
0
0
0
0
Heat from cleaned waste water and district heating network's returning water. Motiva 2014
11
Kellosaari back-up plant
Large fluidized bed
0.3 - 0.5
0
0
0
0
0
-1
0
Only produces electric power
12
Kymijoki River's plants
None
1
0
0
0
0
0
0
0
Hydropower
13
Loviisa nuclear heat
None
0
-0.4 - -0.1
1
0
0
0
0
0
Assumes that for each MWh heat produced, 0.1-0.2 MWh electricity is lost in either production or when heat is pumped to Helsinki.
14
Neste oil refinery heat
None
0
-0.31 - -0.27
1
0
0
0
0
0
Motiva 2014
15
Salmisaari A&B
Large fluidized bed
0.32
0
0.59
0
-1
0
0
0
Capacity: electricity 160 MW heat 300 MW loss 46 MW
16
Sea heat pump
None
0
-0.36 - -0.31
1
0
0
0
0
0
Motiva 2014
17
Sea heat pump for cooling
None
0
-0.36 - -0.31
0
1
0
0
0
0
Assuming the same as for heating
18
Small-scale wood burning
Household
0
0
0.5 - 0.9
0
0
0
0
-1
19
Small gas heat plants
Large fluidized bed
0
0
0.91
0
0
-1
0
0
20
Small fuel oil heat plants
Large fluidized bed
0
0
0.91
0
0
0
-1
0
21
Suvilahti power storage
None
1
0
0
0
0
0
0
0
22
Suvilahti solar
None
1
0
0
0
0
0
0
0
23
Vanhakaupunki museum
None
1
0
0
0
0
0
0
0
Hydropower
24
Vuosaari A
Large fluidized bed
0.455
0
0.455
0
0
-1
0
0
Capacity: electricity 160 MW heat 160 MW loss 30 MW
25
Vuosaari B
Large fluidized bed
0.5
0
0.41
0
0
-1
0
0
Capacity: electricity 500 MW heat 424 MW loss 90 MW
26
Vuosaari C biofuel
Large fluidized bed
0.47
0
0.44
0
0
0
0
-1
27
Wind mills
None
1
0
0
0
0
0
0
0
Notes about the data in the table:
Household air heat pumps data from heat pump comparison[2]
Household geothermal heat data from Energy Department of the United States: Geothermal Heat Pumps[3]
Small-scale wood burning data from Energy Department of the United States: Wood and Pellet Heating[4]
Loss of thermal energy through distribution is around 10 %. From Norwegian Water Resources and Energy Directorate: Energy in Norway.[5]
Sustainable Energy Technology at Work: Use of waste heat from refining industry, Sweden.[6]
Chalmers University of Technology: Towards a Sustainable Oil Refinery, Pre-study for larger co-operation projects[7]
CHP diesel generators are regular diesel generators, but they are located in apartment houses and operated centrally. This way, it is possible to produce electricity when needed and use the excess heat, instead of district heat, to warm up the hot water of the house.
Motiva estimates for heat pumps processes and costs for heating:[8]
Mechanical heat pumps usually have COP (coefficient of performance, thermal output energy per electric input energy needed) is 2.5 - 7.5.
In district heating, mechanical heat pumps have typically COP around 3.
Absorption heat pumps have COP typically 1.5 - 1.8. They do not use much electricity but they need either hot water or steam to operate. Therefore, they are not suitable for producing district heat from warm water with temperatures in the range of 25 - 30 °C (Neste) or 10-15 °C (sea heat).
The report uses these values for energy prices (€/MWh): bought electricity 50, process steam 25, wood chip 20, district heating 40, own excess heat 0.
The investment cost of a heat pump system (ominaiskustannus) in the cases described in this report were 0.47-0.73 M€/MWth for mechanical heat pumps and 0.072 - 0.102 M€/MWth for absorption heat pumps. These values do not include the pipelines needed, which may vary a lot; in these cases the pipeline costs were 0.1 - 2.5 times the cost of the heat pump.
The energy efficiency is theoretically COP = Tout / (Tout - Tin), and the actual COP values are typically 65 - 75 % of that. If we assume that we want 95 °C district heat out, we get
for sea heat pumps: COP = 368 K / (368 K - 283 K) = 4.3 ideally and in practice 2.8 - 3.2. Electricity needed per 1 MWh output: 0.31 - 0.36 MWh.
Neste process heat: COP = 368 K / (368 K - 303 K) = 5.7 ideally and in practice 3.7 - 4.2. Electricity needed per 1 MWh output: 0.23 - 0.27 MWh (plus what is needed for pumping the heat for 25 km, say + 0.04 MWh)
## This code is Op_en7311/energyProcess on page [[Helsinki energy production]]
library(OpasnetUtils)
energyProcess <- Ovariable("energyProcess", ddata = "Op_en7311", subset = "Heat, power and cooling processes")
# Add Time index so we can make Time dependant decisions
energyProcess@data <- merge(energyProcess@data, data.frame(Time = 1880:2080))
objects.store(energyProcess)
cat("Ovariable energyProcess stored.\n")
Plant specifications
These equations below aim to reflect the energy production facilities and capabilities. The min and max values tell about the range of energy production of the plant, and the cost values tell the costs of building and running the powerplant.
Note! Maintenance cost only contains costs that do not depend on activity. Operational cost contains costs that depend on activity but NOT fuel price; those are calculated separately based on energy produced.
biofuels (pellets, wood chips and possibly biochar)
2
2025-2070
CHP diesel generators
0
1441
144
1
1
Assuming all of Helsinki's apartment houses were fitted with 100 kW generators.
3
2025-2080
Deep-drill heat
0
300
300-900
9.6
40
Investment cost from ETSAP
4
1965-2040
Hanasaari
0
640
0
9.6
8
95% coal, 5% pellets. Assume cost of running and maintenance in coal plants 15€/kW (Sähköenergian kustannusrakenne)
5
2010-2060
Household air heat pumps
0
112
200-300
10
5
Assuming all of Helsinki's detached and row houses were fitted with air heat pumps
6
2010-2060
Household air conditioning
0
67
150-200
10
5
7
2016-2060
Household geothermal heat
0
335
380-450
10
5
Assuming all of Helsinki's detached and row houses were fitted with geothermal heat pumps
8
2020-2035
Household solar
0
105
220-250
5
5
Assuming 700000 m2 suitable for solar panels.
9
2010-2070
Katri Vala cooling
0
60
0
10
3
waste water. Max from Helen
10
2005-2065
Katri Vala heat
0
90
0
10
3
waste water. Max from Helen
11
1980-2050
Kellosaari back-up plant
0
120
0
10
20
oil
12
1980-2070
Kymijoki River's plants
0
60
0
10
1-4
hydropower
13
2022-2080
Loviisa nuclear heat
0
1800-2600
400-4000
10
5
Investment cost includes energy tunnel (double of Neste) but NOT building cost of plant. Some estimate for typical district heat pipes on ground is 2 M€/km; this is clearly a minimum for this project.
14
2020-2060
Neste oil refinery heat
0
300
200-500
10
5
15
1975-2050
Salmisaari A&B
0
506
0
7.6
8
95% coal, 5% pellets
16
2020-2070
Sea heat pump
0
225
280
10
4
17
2020-2070
Sea heat pump for cooling
0
225
280
10
4
18
1980-2070
Small-scale wood burning
78
78
0
1
0
Assuming 70% of Helsinki's detached and row houses have a working fireplace. Operation costs for consumer assumed to be 0.
19
1980-2070
Small gas heat plants
0
600
0
5
5
20
1980-2070
Small fuel oil heat plants
0
1600
0
5
5
21
2015-2040
Suvilahti power storage
-1.2
1.2
100
10
5
electricity storage 0.6 MWh
22
2013-2070
Suvilahti solar
0
0.34
0
10
5
23
1880-2070
Vanhakaupunki museum
0
0.2
0
10
0
water
24
1991-2070
Vuosaari A
0
320
0
5
5
natural gas
25
1998-2070
Vuosaari B
0
924
0
5
5
natural gas
26
2018-2070
Vuosaari C biofuel
0
1331
650
10
9
80-100% biofuels, rest coal
27
2017-2060
Wind mills
0
10
12
0.07-0.15
7-13
upper limit from EWEA-report: The economics of wind energy
28
2016-2070
Data center heat
0
150
70.5-109.5
5
0
Investment cost 0.47-0.73 M€/MWth based on Motiva 2014. Cooling is needed anyway, so assumes operation costs to be 0.
### This code is Op_en7311/nondynsupply on page [[Helsinki energy production]].
### nondynsupply is about non-dynamic energy supply. In other words, it cannot be optimised based on need.
library(OpasnetUtils)
# [[Helsinki energy supply]]
nondynsupply <- Ovariable("nondynsupply", ddata = "Op_en7311", subset = "Non-adjustable energy production")
objects.store(nondynsupply)
cat("Ovariable nondynsupply stored.\n")
Fuel availability
Wood
The byproducts of forest industry make up the bulk of fuel wood, and its quantity is almost completely dependent of the production of the forest industry's main products. Therefore it makes sense to calculate the amount of fuel wood usable in the future using the
predictions about the volume of forest industry's production in coming years.
For example, the maximum potential production of woodchips is calibrated so, that it will reach 25 TWh in year 2020, and it is expected slowly increase to 33 TWh by year 2050. The production potential for firewood (for small scale heating) is expected to remain about
the same at just under 60 PJ. The import of wood fuels is estimated to be 3 TWh at most. [21]
Fuel use by heating type
Helsinki-specific data about connections between Heating and fuel usage. Generic data should be taken from Energy balance. Because all Helsinki-specific data is given in the energyProcess table, this only contains dummy data.
## This code is Op_en4151/ on page [[Prices of fuels in heat production]].
library(OpasnetUtils)
objects.latest("Op_en4151", code_name="fuelPrice")
oprint(EvalOutput(fuelPrice))
Rationale
This page contains prices for electricity, district heating, liquid fuels and consumer prices of hard coal, natural gas and domestic fuels in heat production in Finland. all data is based on knowledge of Statistics Finland In the end of data section there will be also data for maintenance and investment costs!!
This code calculates the price of fuels including tax.
library(OpasnetUtils)
fuelPrice <- Ovariable(
"fuelPrice",
dependencies = data.frame(
Name = c("fuelPriceRaw", "fuelTax"),
Ident = c("Op_en4151/fuelPriceRaw", "Op_en4151/fuelTax")
),
formula = function(...){
out <- fuelPriceRaw
out$Result <- 0 #result(out)
out <- merge(out, Ovariable(output = data.frame(Tax = unique(fuelTax$Tax)), marginal = TRUE))
time <- as.POSIXct(strptime(paste("1.1.", as.character(out$Time), sep = ""), format = "%d.%m.%Y"))
for (i in 1:nrow(out@output)) {
temp <- merge(out[i,colnames(out@output) != "Time"], fuelTax)
if (nrow(temp@output) > 0) out$Result[i] <- temp$fuelTaxResult[which.min(abs(temp$Time - time[i]))]
}
out <- oapply(out, NULL, sum, "Tax")
out <- out + fuelPriceRaw
levels(out$Fuel)[levels(out$Fuel) == "Electricity"] <- "Electricity_taxed"
out@output <- orbind(out, fuelPriceRaw[fuelPriceRaw$Fuel == "Electricity",])
result(out)[is.na(result(out))] <- out$fuelPriceRawResult[is.na(result(out))]
return(out)
}
)
objects.store(fuelPrice)
Prices of fuels without tax
A previous version was based on several info sources [8]. However, this resulted in inconsistent prices between fuels. Now the table is simplified and made more robust using IEA estimates VTT-R-03704-14_YDINPAP_(1) page 8, with the cost of accuracy. There should be an expert panel discussing the interdependencies of fuel prices to get the correct order.
Coal, Gas/Natural Gas, Light oil and Bio/Biofuel and Peat data from Statistics Finland[22]. Bio/Biofuel stands for wood chips.
Oil prices from U.S. Energy Information Administration[23]. All values converted first to 2015 dollars (inflation correction) and then to 2015 euros using the current (21.7.2015) exchange ratio of 1 $ = 0.923 €. Price per barrel then converted to price per MWh presuming that the energy released by burning one barrel is about 5.8 x 106 BTU = 1.7 MWh, U.S. Internal Revenue Service[24].
Compound Annual Growth Rates (CAGR) for Coal, Gas/Natural Gas, Oil, Crude oil, Fuel oil, Heavy oil and Light oil from U.S. Energy Information Administration's "Annual Energy Outlook 2015 with projections to 2040" [25]. CAGR for Bio/Biofuel and Peat estimated to be 2 %, based on the price history.
District heating data is based on the price for an apartment building (volume 10 000 m3, energy need 450 MWh/a). Data from Statistics Finland [10]. We assume that the price of district cooling is the same (although district cooling has been available only for a few years).
Uncertainties are assumed in the same way for all fuels: a certain percentage up and down from the best estimate, varying by year. The following uncertainties were used: 2025 5 %, 2035 10 %, 2045 15 %, 2055 20 %, 2065 30 %. It should be noted that uncertainties are different for different fuels, but we were unable to estimate fuel-specific uncertainties.
## This is code Op_en4151/fuelPriceRaw on page Prices of fuels in heat production
library(OpasnetUtils)
fuelPriceRaw <- Ovariable("fuelPriceRaw", ddata = "Op_en4151", subset = "Fuel prices")
objects.store(fuelPriceRaw)
cat("Ovariable fuelPriceRaw stored.\n")
Summing up the amount of energy produced and subtracting the amount of energy consumed within a time period gives the energy balance. Since the electricity grid and district heat network lack significant storage mechanics, the balance has to be virtually zero over short periods. When considering the balance of a particular area (e.g. Helsinki), we can make the assumption that electricity can be imported and exported in international markets. The energy in the district heat network, however, has to be produced locally. This sets up the non-trivial problem of optimising production so that there are no significant deficits as well as minimising losses and maximising profits. This problem is solved (to some extent) by market forces in the real world.
In Opasnet, there are two different ways to calculate energy balance. Our most recent energy balance model uses linear programming tools to solve an optimum for the activity of a given set of production units in simulated instances created by the main model. The main model is responsible for the decision making aspects, while the energy balance optimisation only functions as an approximation of real world market mechanics. This version was used in Helsinki energy decision 2015.
The previous version was based on setting up a set of linear equations describing the inputs, outputs, and shares of different energy and plant processes. This approach is less flexible, because it does not use an optimising function and everything must be described as linear (or piecewise linear). However, this approach was successfully used in Energy balance in Kuopio and Energy balance in Suzhou.
Example of energy balance model: District heat flow in Helsinki. The scenarios are from the assessment Helsinki energy decision 2015.
Example of energy balance model: Electric power heat flow in Helsinki. The scenarios are from the assessment Helsinki energy decision 2015.
Example of energy balance model: Incomes and costs of energy production in Helsinki. The scenarios are from the assessment Helsinki energy decision 2015.
This code is an example how the energy balance model is used in a city case. The data comes from Helsinki energy decision 2015.
library(OpasnetUtils)
library(ggplot2)
openv.setN(0) # use medians instead of whole sampled distributions
# Download case-specific data, in this case from Helsinki.
objects.latest("Op_en7237", code_name = "intermediates") # [[Helsinki energy decision 2015]]
# Download energy balance model and its parts:
# EnergyConsumerDemandTotal
# EnergyFlowHeat
# EnergyFlowOther
# EnergyNetworkDemand
# EnergyNetworkOptim
# fuelUse
# EnergyNetworkCost
objects.latest("Op_en5141", code_name = "EnergyNetworkOptim") # [[Energy balance]]
oprint(summary(EvalOutput(EnergyNetworkOptim)))
Rationale
Energy balance with linear programming
The linear programming problem is set up as follows.
For each production unit: let xi be activity of the plant. Lets also have variables yj for deficits and excesses for each type of energy produced.
The objective function is the function we are optimising. Each production unit has a unit profit per activity denoted by ai which is determined by the amount of different input commodities (e.g. coal) per amount of different output commodities (i.e. electricity and heat) and their market prices. Also, lets say we want to make sure that district heat demand is always met when possible and have a large penalty factor for each unit of heat demand not met (1 M€ in the model). In addition, it must be noted that excess district heat becomes wasted so it counts as loss. Let these deficit and excess related losses be denoted by bj. The whole objective function then becomes: sum(xiai) + sum(yjbj).
The values of variables are constrained by equalities and inequalities: the sum of production of a commodity is equal to its demand minus deficit plus excess, activity is constrained by the maximum capacity and all variables are non-negative by definition.
This can be efficiently solved by computers for each given instance. Production wind-up and wind-down is ignored, since time continuity is not considered. As a consequence fuel limits (e.g. diminishing hydropower capacity) are not modelled completely either.
## This code is Op_en5141/EnergyNetworkOptim [[Energy balance]].
library(OpasnetUtils)
EnergyConsumerDemandTotal <- Ovariable("EnergyConsumerDemandTotal",
dependencies = data.frame(
Name = c(
"EnergyConsumerDemand"
)
),
formula = function(...) {
energy <- oapply(EnergyConsumerDemand, NULL, sum, c("Renovation", "Efficiency"))
energy <- energy * 1e-6 # W to MW
return(energy)
}
)
EnergyFlowHeat <- Ovariable("EnergyFlowHeat",
dependencies = data.frame(
Name = c(
"EnergyConsumerDemandTotal",
"fuelShares"
)
),
formula = function(...) {
EnergyFlowHeat <- EnergyConsumerDemandTotal[EnergyConsumerDemandTotal$Consumable %in% c("Heating", "Hot water"),] * fuelShares
return(EnergyFlowHeat)
}
)
EnergyFlowOther <- Ovariable("EnergyFlowOther",
dependencies = data.frame(
Name = c(
"EnergyConsumerDemandTotal",
"temperene",
"nontemperene"
)
),
formula = function(...) {
EnergyFlowOther <- merge(
EnergyConsumerDemandTotal[!EnergyConsumerDemandTotal$Consumable %in% c("Heating", "Hot water"),],
unique(OpasnetUtils::combine(temperene, nontemperene)[,c("Consumable", "Fuel")])
)
EnergyFlowOther@name <- EnergyConsumerDemandTotal@name
return(EnergyFlowOther)
}
)
EnergyNetworkDemand <- Ovariable("EnergyNetworkDemand",
dependencies = data.frame(
Name = c(
"EnergyFlowHeat",
"EnergyFlowOther",
"EnergyConsumerDemandTotal"
)
),
formula = function(...) {
demand <- OpasnetUtils::combine(
EnergyFlowHeat[EnergyFlowHeat$Fuel %in% c("Heat", "Electricity"), colnames(EnergyFlowHeat@output) != "Burner"],
EnergyFlowOther[EnergyFlowOther$Fuel %in% c("Cooling", "Electricity"),]
)
# EnergyFlowOther should not have other flows
#assign("EnergyFlowOther", NULL, .GlobalEnv)
#assign("EnergyFlowHeat", EnergyFlowHeat[!EnergyFlowHeat$Fuel %in% c("Heat", "Electricity"),], .GlobalEnv)
# Total consumer demand also unnecessary?
#assign("EnergyConsumerDemandTotal", NULL, .GlobalEnv)
NAindex <- sapply(lapply(lapply(demand@output, unique), is.na), sum)
NAindex <- names(NAindex)[NAindex != 0]
demand@output <- fillna(demand@output , NAindex)
demand <- oapply(demand, NULL, sum, c("Heating", "Consumable"))#, "Renovation", "Efficiency")) # Assumes no NA
return(demand)
}
)
EnergyNetworkOptim <- Ovariable("EnergyNetworkOptim",
dependencies = data.frame(
Name = c(
#"energy", # Energy supply and demand in the system
#"nondynsupply", # Energy supply that is not optimised.
#"requiredName", # Fuel name for the energy type that is balanced for inputs and outputs.
#"fuelShares", # Shares of fuels for different heating types.
"energyProcess", # Matrix showing inputs and outputs of the process of each plant.
"fuelPrice", # Prices of fuels (by Time)
#"timelylimit", # Maximum energy supply in the given conditions at a particular time point
"plantParameters", # Capacities and costs of running each plant
#"temperene",
#"nontemperene",
"EnergyNetworkDemand"
),
Ident = c(
#NA,
#NA,
#NA,
#NA,
NA,
"Op_en4151/fuelPrice",
#NA,
NA,
#"Op_en5488/temperene",
#"Op_en5488/nontemperene",
NA
)
),
formula = function(...) {
require(lpSolve)
require(plyr)
require(reshape2)
# Prices are given as €/MWh while plant activity is given as MW while optimisation is done
# with time resolution = 1 day
optim_time_period_adj <- 24
# List non-interating indices
exclude <- c("Fuel", "Plant", "Parameter", "Heating", "Consumable", "Burner")
# Calculate unitprofit for plant activity energy flows
unitp <- fuelPrice * energyProcess * optim_time_period_adj
unitp <- oapply(unitp, NULL, sum, c("Fuel", "Burner"))
colnames(unitp@output)[colnames(unitp@output) == "Result"] <- "UnitPrice"
# Reshape crucial variables to reduce merging difficulty
ePmargs <- colnames(energyProcess@output)[energyProcess@marginal]
energyProcess@output <- dcast(
energyProcess@output,
as.formula(paste(paste(ePmargs[ePmargs != "Fuel"], collapse = "+"), "~ Fuel")),
value.var = paste(energyProcess@name, "Result", sep = ""),
fill = NA
)
energyProcess@marginal <- colnames(energyProcess@output) %in% ePmargs
pPmargs <- colnames(plantParameters@output)[plantParameters@marginal]
plantParameters@output <- dcast(
plantParameters@output,
as.formula(paste(paste(pPmargs[pPmargs != "Parameter"], collapse = "+"), "~ Parameter")),
value.var = paste(plantParameters@name, "Result", sep = ""),
fill = NA
)
plantParameters@marginal <- colnames(plantParameters@output) %in% pPmargs
EnergyNetworkDemand$Fuel <- paste(EnergyNetworkDemand$Fuel, "Demand", sep = "")
ENDmargs <- colnames(EnergyNetworkDemand@output)[EnergyNetworkDemand@marginal]
EnergyNetworkDemand@output <- dcast(
EnergyNetworkDemand@output,
as.formula(paste(paste(ENDmargs[ENDmargs != "Fuel"], collapse = "+"), "~ Fuel")),
value.var = paste(EnergyNetworkDemand@name, "Result", sep = ""),
fill = NA
)
EnergyNetworkDemand@marginal <- colnames(EnergyNetworkDemand@output) %in% ENDmargs
fuelPrice$Fuel <- paste(fuelPrice$Fuel, "Price", sep = "")
fPmargs <- colnames(fuelPrice@output)[fuelPrice@marginal]
fuelPrice@output <- dcast(
fuelPrice@output,
as.formula(paste(paste(fPmargs[fPmargs != "Fuel"], collapse = "+"), "~ Fuel")),
value.var = paste(fuelPrice@name, "Result", sep = ""),
fill = 0
)
fuelPrice@marginal <- colnames(fuelPrice@output) %in% fPmargs
# Duplicates a couple of values, but provides better performance
vars <- merge(EnergyNetworkDemand, plantParameters)
vars <- merge(vars, energyProcess)
vars <- merge(vars, fuelPrice)
vars <- merge(vars, unitp)
# Split into iterations
#vars <- dlply(vars@output, colnames(vars@output)[vars@marginal & !colnames(vars@output) %in% exclude], I)
#out <- data.frame()
include <- union(ePmargs, pPmargs)
include <- union(include, ENDmargs)
include <- union(include, fPmargs)
include <- include[!include %in% exclude]
# Optimisation iteration function
optf <- function(varsi) {
# Acquire relevant sections of variables with respect to index iteration
demandi <- unlist(varsi[1,grepl("Demand$", colnames(varsi))])
names(demandi) <- gsub("Demand$", "", names(demandi))
fuelPricei <- unlist(varsi[1,grepl("Price$", colnames(varsi))])
names(fuelPricei) <- gsub("Price$", "", names(fuelPricei))
fuelPricei <- fuelPricei[names(fuelPricei) != "Unit"]
plants <- as.character(unique(varsi[["Plant"]])) # Used for ordering parameters
if (nrow(varsi) != length(plants)) stop("Wrong subset given to energy production optimizer! Each row not unique with respect to Plant.")
lower <- varsi[["Min"]]
upper <- varsi[["Max"]]
opcost <- varsi[["Operation cost"]]
unitpi <- varsi[["UnitPrice"]]
#commodities <- as.character(unique(energyProcessi$Fuel)) # Used for ordering parameters
# Constraint matrix
# Total number of variables:
# * plant activity for each plant
# * (commodity flow total)
# * commodity deficit and excess with respect to demand
# * opeartion costs
nvars <- length(plants) + length(demandi) * 2 + 1
# Rows:
# * demand rows and their constraints
# * (commodity flow rows)
# * plant activity constraints
# * operation costs row
nrows <- 3 * length(demandi) + 2 * length(plants) + 1
M <- matrix(
0,
nrows,
nvars
)
sign <- rep(">=", nrows)
constant <- rep(0, nrows)
for (j in 1:length(demandi)) {
# Total flow (production - consumption between all processes) + deficit = demand + excess
# which translates to:
# total_commodity_flow + deficit - excess == constant ("static" demand)
M[j, length(plants) + j*2 + 1:2 - 2] <- c(1, -1)
sign[j] <- "=="
constant[j] <- demandi[j]
# Plant flows by activity
M[j, 1:length(plants)] <- varsi[[names(demandi)[j]]]
}
rowi <- length(demandi)
diag(M[
(rowi + 1):(rowi + length(plants)),
1:length(plants)
]) <- 1
diag(M[
(rowi + length(plants) + 1):(rowi + 2 * length(plants)),
1:length(plants)
]) <- 1
constant[rowi + 1:length(plants)] <- lower
constant[rowi + length(plants) + 1:length(plants)] <- upper
sign[rowi + length(plants) + 1:length(plants)] <- "<="
rowi <- rowi + 2 * length(plants)
# Excess and deficit must both be positive
diag(M[
rowi + 1:(2 * length(demandi)),
length(plants) + #length(commodities) +
1:(2 * length(demandi))
]) <- 1
# Operating cost
M[nrows, 1:length(plants)] <- opcost * optim_time_period_adj
M[nrows, nvars] <- -1
sign[nrows] <- "=="
# Objective function (profits - penalty)
objective <- rep(0, nvars)
# Commodity prices
#objective[length(plants) + (1:length(commodities))[!is.na(fuelPricej)]] <- result(fuelPricei)[fuelPricej[!is.na(fuelPricej)]]
# Penalize heat deficit heavily,
# electricity deficit has no impact on profit assuming company sells to customers at the purchase price
# cooling deficit has no effect
if ("Heat" %in% names(demandi)) {
objective[length(plants) + match("Heat", names(demandi)) * 2 - 1] <- -1e6
}
# Electricity excess can be sold in the markets,
# heat and cooling (and other?) excess, however, cannot, hence they need to be deducted from profits
objective[
length(plants) + which(names(demandi) != "Electricity") * 2
] <- - fuelPricei[match(names(demandi)[names(demandi) != "Electricity"], names(fuelPricei))] * optim_time_period_adj
# Operation cost
objective[nvars] <- -1
# Unit profit
objective[1:length(plants)] <- unitpi
# Do actual linear programming
temp <- lp(
"max",
objective,
M,
sign,
constant
)
# Format results
outi <- varsi[rep(1, nvars + 1), include, drop = FALSE]
outi[["Process_variable_type"]] <- rep(
c("Activity", "Flow", "Misc"),
c(length(plants), nvars - length(plants) - 1, 2)
)
outi[["Process_variable_name"]] <- c(
as.character(plants),
paste(rep(names(demandi), each = 2), c("Deficit", "Excess"), sep = ""),
"Operation cost",
"Profit"
)
outi[["Result"]] <- c(temp$solution, temp$objval)
return(outi)
}
#aggregate(out$Result[out$Process_variable_name == "Profit"], out[out$Process_variable_name == "Profit",colnames(out) %in% c("Time", "Temperature") | grepl("Policy", colnames(out))], mean)
#aggregate(out$Result[out$Process_variable_name == "HeatDeficit"], out[out$Process_variable_name == "HeatDeficit", colnames(out) %in% c("Time", "Temperature") | grepl("Policy", colnames(out))], mean)
#aggregate(out$Result[out$Process_variable_name == "Hanasaari"], out[out$Process_variable_name == "Hanasaari", colnames(out) %in% c("Time", "Temperature") | grepl("Policy", colnames(out))], mean)
# Error handling
optfsecure <- function(varsi) {
ret <- tryCatch(
optf(varsi),
error = function(e) return(NULL)
)
if (!is.null(ret)) {
return(ret)
} else {
warning(paste("EnergyNetworkOptim failed optimising a permutation with error:", geterrmessage()))
return(data.frame())
}
}
out <- ddply(vars@output, colnames(vars@output)[vars@marginal & !colnames(vars@output) %in% exclude], optfsecure)
out <- Ovariable(output = out, marginal = !grepl("Result$", colnames(out)))
return(out)
}
)
fuelUse <- Ovariable("fuelUse",
dependencies = data.frame(
Name = c(
"EnergyNetworkOptim",
"EnergyNetworkDemand",
"energyProcess",
"EnergyFlowHeat",
"temperdays"
)
),
formula = function(...){
# Calculate flows resulting from plant activity
EnergyFlow <- EnergyNetworkOptim[EnergyNetworkOptim$Process_variable_type == "Activity", colnames(EnergyNetworkOptim@output) != "Process_variable_type"]
colnames(EnergyFlow@output)[colnames(EnergyFlow@output) == "Process_variable_name"] <- "Plant"
EnergyFlow <- EnergyFlow * energyProcess
# Realised consumption of centrally produced commodities (demand - deficit or prod - excess)
# NOTE: centralized flows which consume electricity for example should be separated
# from consumer demand and consumption, so prod - excess is not accurate
deficit <- EnergyNetworkOptim[
grepl("Deficit$", EnergyNetworkOptim$Process_variable_name),
colnames(EnergyNetworkOptim@output) != "Process_variable_type"
]
colnames(deficit@output)[colnames(deficit@output) == "Process_variable_name"] <- "Fuel"
deficit$Fuel <- gsub("Deficit$", "", deficit$Fuel)
# Assume that the power company purchases any electricity deficit from the markets in order to satisfy demand
result(deficit)[deficit$Fuel == "Electricity"] <- 0
#excess <- EnergyNetworkOptim[
# grepl("Excess$", EnergyNetworkOptim$Process_variable_name),
# colnames(EnergyNetworkOptim@output) != "Process_variable_type"
#]
#colnames(excess@output)[colnames(excess@output) == "Process_variable_name"] <- "Fuel"
#excess$Fuel <- gsub("Excess$", "", excess$Fuel)
#temp <- oapply(EnergyFlow, NULL, sum, c("Plant", "Burner"))
#real <- oapply(EnergyFlow[EnergyFlow$Fuel %in% unique(excess$Fuel)], NULL, sum, c("Plant", "Burner"))
#real <- real - excess
real <- EnergyNetworkDemand - deficit
# Combine
real$Burner <- "None"
real$Plant <- "Domestic"
heating <- EnergyFlowHeat[!EnergyFlowHeat$Fuel %in% unique(deficit$Fuel),]
heating <- oapply(heating, NULL, sum, "Consumable")
heating$Plant <- "Domestic"
heating$Heating <- NULL
EnergyFlow <- OpasnetUtils::combine(0 - EnergyFlow, real, heating)
EnergyFlow <- unkeep(EnergyFlow, sources = TRUE, prevresults = TRUE)
#EnergyFlowTest <- oapply(EnergyFlow, c("Time", "Temperature", "Fuel"), sum)
#ggplot(EnergyFlowTest@output, aes(x = Temperature, y = Result, group = Fuel, color = Fuel)) + geom_line() + facet_wrap( ~ Time)
EnergyFlow <- EnergyFlow * temperdays * 3600 * 24
EnergyFlow <- oapply(EnergyFlow, cols = c("Temperature"), FUN = sum)
return(EnergyFlow)
}
)
EnergyNetworkCost <- Ovariable("EnergyNetworkCost",
dependencies =
data.frame(
Name = c(
"plantParameters",
"EnergyNetworkOptim",
"temperdays"
),
Ident = c(
NA,
"Op_en5141/EnergyNetworkOptim", # [[Energy balance]]
NA
)
),
formula = function(...) {
oper <- plantParameters[plantParameters@output$Parameter == "Max" , colnames(plantParameters@output) != "Parameter"]
result(oper)[result(oper) != 0] <- 1
oper <- plantParameters * oper
# Take the first year when a plant is operated and put all investment cost there.
investment <- oper[oper@output$Parameter == "Investment cost" , colnames(oper@output) != "Parameter"]
investment <- investment[result(investment) > 0 , ]
investment <- investment[order(investment@output$Time) , ]
investment <- investment[!duplicated(investment@output[investment@marginal & colnames(investment@output) != "Time"]) , ]
investment <- unkeep(investment, sources = TRUE)
#investment <- oapply(investment, cols = "Plant", FUN = sum)
maintenance <- oper[oper@output$Parameter == "Management cost" , colnames(oper@output) != "Parameter"]
maintenance <- unkeep(maintenance, sources = TRUE)
#maintenance <- oapply(maintenance, cols = c("Plant"), FUN = sum)
operation <- EnergyNetworkOptim[EnergyNetworkOptim@output$Process_variable_name == "Operation cost" , ]
operation <- operation * temperdays * 10 * 1E-6 # For 10-year periods, € -> M€
operation <- oapply(operation, cols = c("Temperature"), FUN = sum)
operation <- unkeep(operation, cols = c("Process_variable_name", "Process_variable_type"), sources = TRUE, prevresults = TRUE)
operation <- operation * Ovariable(output = data.frame(Plant = "Operation", Result = 1), marginal = c(TRUE, FALSE))
cost <- OpasnetUtils::combine(investment, maintenance, operation)
marginals <- character()
for(i in colnames(cost@output)[cost@marginal]) {
if(any(is.na(cost@output[[i]]))) marginals <- c(marginals, i)
}
if(length(marginals) > 0) {
cost@output <- fillna(cost@output, marginals)
warning(paste("In combine had to fillna marginals", marginals, "\n"))
}
return(cost)
}
)
objects.store(EnergyConsumerDemandTotal, EnergyFlowHeat, EnergyFlowOther, EnergyNetworkDemand, EnergyNetworkOptim, fuelUse, EnergyNetworkCost)
cat("Ovariables EnergyConsumerDemandTotal, EnergyFlowHeat, EnergyFlowOther, EnergyNetworkDemand, EnergyNetworkOptim, EnergyNetworkCost and fuelUse stored.\n")
Fuel use and fuel shares in generic processes
There is an alternative way for calculating fuel use. It is based on the idea that for each heating type, there is a constant share of fuels used. For some heating types, this is generic and is shown on this page. For some others, the constant is case-specific and is determined on a case-specific page.
The table below contains connections of heating types and fuel usage in generic situations. There may be case-specific differences, which must be handled separately.
## This code is Op_en2719/ on page [[Emission factors for burning processes]].
library(OpasnetUtils)
library(ggplot2)
objects.latest("Op_en2791", code_name = "emissionstest")
objects.latest("Op_en2791", code_name = "emissionFactors")
oprint(summary(EvalOutput(emissionFactors)))
Rationale
HSY emission factors
This section describes the use of HSY emission factors used for CO2 emission estimates. This should be merged with the other data, but it is first described as such.
These emission factors are derived from the Helsinki Region Environmental Services (HSY) climate data [11] on 27th Nov 2018. Emission factors were initially calculated for every entry in the data file (4180 rows in total). Then, emission sectors were compared along timeline for each city separately. We noticed that several sectors shared practically the same emission factors with few differences, which are not plausible as real differences (Table Sector classification). So, we took a mean of the values to represent all sectors in the same EFclass group.
Then, we compared results from different cities. There were minor changes that may be due to some real differences in data, and also a some changes that looked like artefact. However, there were two sectors, namely district heating and fuels that clearly differed between cities in a plausible way. Therefore, we used city-specific emission factors for these two sectors and those of Helsinki for all other sectors (because Helsinki seemed to have the least amount of artefact in the data).
# This is code Op_en2791/HSYdata on page [[Emission factors for burning processes]]
library(OpasnetUtils)
HSYdata <- Ovariable("HSYdata",ddata="Op_en2791",subset="HSY CO2 emission and energy consumption")
EFclass <- Ovariable("EFclass", ddata="Op_en2791", subset="Sector classification")
objects.store(HSYdata, EFclass)
cat("Ovariables HSYdata, EFclass stored.\n")
###This code is Op_en2791/emissionstest on page [[Emission factors for burning processes]].
library(OpasnetUtils)
####### Calculate emissions
emissions <- Ovariable("emissions",
dependencies = data.frame(
Name = c(
"fuelUse", # Use of fuels (in MJ)
"emissionFactors", # Emissions of different pollutants per energy (mg /MJ)
"emissionLocations" # Locations where the emissions are emitted, indexed by City_area or Plant
),
Ident = c(
NA, # typically from [[Energy balance]] but there are other ways to calculate
"Op_en2791/emissionFactors", # [[Emission factors for burning processes]]
NA # Case-specific data on locations of emissions.
)
),
formula = function(...) {
# convert from mg to ton emitted.
out <- fuelUse * emissionFactors * 1E-9
if(!"City_area" %in% colnames(out@output)) {
out <- out * Ovariable(
output = data.frame(City_area = "Not known", Result = 1),
marginal = c(TRUE, FALSE)
)
}
out <- unkeep(out * emissionLocations, sources = TRUE, prevresults = TRUE)
out@output$Emission_site <- as.factor(ifelse(
out@output$Emission_site == "At site of consumption",
as.character(out@output$City_area),
as.character(out@output$Emission_site)
))
out <- oapply(
out,
cols = c("Burner", "City_area"),
FUN = sum
)
return(out)
}
)
objects.store(emissions)
cat("Ovariable emissions stored.\n")
Light oil <5 MW Emission factors for burning processes. Light oil 267 kg /MWh
4
Household
Oil
0-10
74200-87222
74200
87222
Light oil <5 MW Emission factors for burning processes. Light oil 267 kg /MWh
5
Household
Other sources
0-10
74200
74200
74200
Same as oil.
6
Household
Coal
0-10
74200-87222
74200
87222
7
Household
Geothermal
0-10
74200-87222
74200
87222
8
Household
Gas
0-3
55650
55650
55650
For PM2.5: one third of that of oil. For CO2: 3/4 of that of oil.
9
Household
Fuel oil
0-10
74200-87222
74200
87222
Light oil <5 MW Emission factors for burning processes. Light oil 267 kg /MWh
10
Domestic
Wood
140 (65.8-263)
74200
0
8333
Other stoves and ovens. Karvosenoja et al. 2008 Just repeat the previous rows to match different wording of burners.
11
Domestic
Biofuel
140 (65.8-263)
74200
0
8333
Other stoves and ovens. Karvosenoja et al. 2008
12
Domestic
Light oil
0-10
74200-87222
74200
87222
Light oil <5 MW Emission factors for burning processes. Light oil 267 kg /MWh
13
Domestic
Oil
0-10
74200-87222
74200
87222
Light oil <5 MW Emission factors for burning processes. Light oil 267 kg /MWh
14
Domestic
Other sources
0-10
74200
74200
74200
Same as oil.
15
Domestic
Coal
0-10
74200-87222
74200
87222
16
Domestic
Geothermal
0-10
74200-87222
74200
87222
17
Domestic
Gas
0-3
55650
55650
55650
For PM2.5: one third of that of oil. For CO2: 3/4 of that of oil.
18
Domestic
Fuel oil
0-10
74200-87222
74200
87222
Light oil <5 MW Emission factors for burning processes. Light oil 267 kg /MWh
19
Diesel engine
Fuel oil
0-10
74200-87222
74200
87222
Light oil <5 MW Emission factors for burning processes. Light oil 267 kg /MWh
20
Diesel engine
Light oil
0-10
74200-87222
74200
87222
21
Diesel engine
Biofuel
0-10
74200-87222
74200
87222
22
Large fluidized bed
Gas
0-3
55650
55650
55650
For PM2.5: one third of that of oil. For CO2: 3/4 of that of oil.
23
Large fluidized bed
Coal
2-20
106000
106000
106000
Same as peat.
24
Large fluidized bed
Wood
2-20
74200
0
74200
Leijupoltto 100-300 MW Emission factors for burning processes. Karvosenoja et al., 2008
25
Large fluidized bed
Biofuel
2-20
74200
0
74200
Leijupoltto 100-300 MW Emission factors for burning processes. Karvosenoja et al., 2008
26
Large fluidized bed
Waste
2-20
74200
0
-50000
CO2trade same as wood. CO2eq is guesswork but it is negative because without burning it would produce methane in landfill
27
Large fluidized bed
Peat
2-20
106000
106000
107500
Leijupoltto 100-300 MW Emission factors for burning processes. Peat 382 kg /MWh
28
Large fluidized bed
Heavy oil
8-22
91111-106000
106000
91111
Leijupoltto 100-300 MW Emission factors for burning processes. Peat 382 kg /MWh
29
Large fluidized bed
Fuel oil
8-22
91111-106000
106000
91111
Leijupoltto 100-300 MW Emission factors for burning processes. Peat 382 kg /MWh
30
Grid
Electricity
1-10
53000
212000
53000
50 % of large-scale burning (because of nuclear and hydro). Heavy oil 279 kg /MWh. Officially, electricity is not CHP but requires a double amount of coal to produce it.
31
None
Electricity_taxed
1-10
53000
212000
53000
50 % of large-scale burning (because of nuclear and hydro). Heavy oil 279 kg /MWh. Officially, electricity is not CHP but requires a double amount of coal to produce it. These emissions are assumed when power plants buy electricity from the grid.
32
None
Electricity
0
0
0
0
We might want to keep these locations in the model, but we assume that emissions are zero.
33
None
Heat
0
0
0
0
We might want to keep these locations in the model, but we assume that emissions are zero.
34
None
Cooling
0
0
0
0
We might want to keep these locations in the model, but we assume that emissions are zero.
Large fluidized bed (Peat) CO2-eq value from Väisänen, Sanni: Greenhouse gas emissions from peat and biomass-derived fuels, electricity and heat — Estimation of various production chains by using LCA methodology[27]
Other CO2-eq values from EKOREM: Sähkölämmitys ja lämpöpumput sähkönkäyttäjinä ja päästöjen aiheuttajina Suomessa.
Classes of climate emissions:
CO2direct
Direct CO2 emissions from the stack
CO2trade
CO2 emissions as they are defined in the emission trade. Non-trade sectors have emission 0.
CO2eq
CO2 emissions as equivalents (i.e. includes methane, N2O and other climate emissions based on life cycle impacts.
In Finland there are about 700 kettles that have under 5MW fuel power. Same amount is between 5 to 50 MW kettles and over 50 MW kettles there are 200 in Finland. One heating power plant can have several kettles. Many 5-50 MW power plants have also less than 5 MW a kettle. [28]
See further discussions about emission factors of wood burning and other topics on the discussion page.D↷
# This is code Op_en2791/emissionFactors on page [[Emission factors for burning processes]].
library(OpasnetUtils)
emissionFactors <- Ovariable(
name = 'emissionFactors',
ddata = 'Op_en2791', # [[Emission factors for burning processes]]
subset = 'Emission factors of energy production'
)
objects.store(emissionFactors)
cat("Objects emissionFactors initiated!\n")
Intake fractions of fine particles
Question
How to calculate exposure based on intake fractions of airborne particulate matter for different emission sources and locations?
Answer
Intake fraction (iF) is the fraction of an emission that is ultimately breathed by someone in the target population. With fine particles, it is often in the range of one in a million, but variation is large. It can be used as a shortcut for calculating exposures in a situation where actual atmospheric fate and transport modelling is not feasible. For fine particles, there is fairly good understanding of the magnitudes of intake fractions in different situations.
[29]
Therefore, they have been successfully used in many assessments.
Intake fraction is defined as
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle iF = \frac{c * P * BR}{E},}
where
iF = intake fraction (unitless after proper unit conversions)
c = exposure concentratíon of the population (µg/m3)
P = population size
BR = breating rate, usually a nominal value 20 m3/d is used
E = emission of fine particles (g/s)
In an assessment, exposure concentration c is solved from the equation and used as exposure in health impact modelling.
## This is code Op_en5813/iFHumbert on page [[Intake fractions of PM]].
library(OpasnetUtils)
iF <- Ovariable("iF", ddata = "Op_en5813", subset = "Intake fractions of PM")
# [[Intake fractions of PM]] Humbert et al 2011 data
colnames(iF@data) <- gsub("[ \\.]", "_", colnames(iF@data))
iF@data$iFResult <- iF@data$iFResult * 1E-6
objects.store(iF)
cat("Ovariable iF (Humbert et al 2011) stored.\n")
Exposure-response functions of environmental pollutants
Question
What are the exposure-response functions (ERF) of environmental pollutants that do not have own pages in Opasnet?
Answer
Exposure-response function (ERF) is a mathematical quantitative description of the relationship between an exposure to an agent and the health responses it causes in the human body. How to actually estimate the responses based on ERF is described in detail on page Health impact assessment.
Relevant example results can be found from here (in Finnish).
## This is code Op_en5827/ on page [[ERFs of environmental pollutants]]
library(OpasnetUtils)
objects.latest("Op_en5827", "initiate")
cat("Exposure-response functions of environmental pollutants.\n")
oprint(summary(EvalOutput(ERF_env)))
Lanphear et al 2005 https://doi.org/10.1289/ehp.7688 using the first increment from 24 to 100 ug/l. Assumes threshold at 24 ug/l although this is strong assumption
# This is code Op_en5827/ERF_env2 on page [[ERF of environmental pollutants]]
# Note! This version has ERF and threshold in the same ovariable.
library(OpasnetUtils)
ERF_env <- Ovariable("ERF_env", ddata = "Op_en5827")
colnames(ERF_env@data) <- gsub(" ", "_", colnames(ERF_env@data))
objects.store(ERF_env)
cat("Ovariable ERF_env stored.\n")
External costs are costs that are not included in the price of a product but still cause negative (or sometimes positive) impacts to the society or stakeholders. The market theory says that if external costs cannot be included in prices e.g. using taxation, the market process will lead to outcomes that deviate from the societal optimum.
Question
What are important external costs in environmental health?
Answer
Health impacts and climate impacts are often not considered in pricing, so they remain as external costs. The code below is used to fetch the variable for models.
## This is code Op_en7379/ on page [[External cost]]
library(OpasnetUtils)
objects.latest("Op_en7379", code_name = "externalities") # [[External cost]] ovariables DALYprice, co2price
cat("Price of a DALY (€)\n")
oprint(summary(EvalOutput(DALYprice)))
cat("Price of a ton of CO2 (€)\n")
oprint(summary(EvalOutput(co2price)))
↑Berntsson T, Persson Elmeroth L, Algehed J, Hektor E, Franck PÅ, Åsblad A, Johnsson F, Lyngfelt A, Gevert B, Richards T: Towards a Sustainable Oil Refinery - Pre-study for larger co-operation projects. Chalmers Energy Centre (CEC) Report 2008:1. Chalmers University of Technology. http://publications.lib.chalmers.se/records/fulltext/69752.pdf
↑ 29.029.1Sebastien Humbert, Julian D. Marshall, Shanna Shaked, Joseph V. Spadaro, Yurika Nishioka, Philipp Preiss, Thomas E. McKone, Arpad Horvath, and Olivier Jolliet. Intake Fraction for Particulate Matter: Recommendations for Life Cycle Impact Assessment (2011). Environmental Science and Technology, 45, 4808-4816.