Energy balance: Difference between revisions

Latest revision as of 07:49, 10 April 2019

This page is a knowledge crystal of subtype method. The page identifier is Op_en5141
Moderator:Jouni (see all)
This page is a stub. You may improve it into a full page.
Upload data {{#opasnet_base_link:Op_en5141}}

Question

What is energy balance and how is it modelled?

Answer

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.

+ Show code - Hide code

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 x_i be activity of the plant. Lets also have variables y_j 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 a_i 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 b_j. The whole objective function then becomes: sum(x_ia_i) + sum(y_jb_j).

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.

Ovariables like EnergyNetworkOptim below are used in Helsinki energy decision 2015. Prices of fuels in heat production are used as direct inputs in the optimising.

+ Show code - Hide code

## 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.

Fuel use in different heating types(-)
Obs	Heating	Burner	Fuel	Fraction	Description
1	Wood	Household	Wood	1
2	Oil	Household	Light oil	1
3	Gas	Household	Gas	1
4	Heating oil	Household	Light oil	1
5	Coal	Household	Coal	1
6	Other sources	Household	Other sources	1
7	No energy source	Household	Other sources	1
8	Geothermal	Grid	Electricity	0.3	Geothermal does not sum up to 1 because more heat is produced than electricity consumed.
9	Centrifuge, hydro-extractor	Grid	Electricity	0.3	Not quite clear what this is but presumably a heat pump.
10	Solar heater/ collector	Grid	Electricity	0.1	Use only; life-cycle impacts omitted.
11	Electricity	Grid	Electricity	1
12	District	Undefined	Heat	1

+ Show code - Hide code


# This is code Op_en5141/fuelSharesgeneric (only generic) on page [[Energy balance]].
library(OpasnetUtils)

fuelSharesgeneric <- Ovariable("fuelSharesgeneric", ddata = "Op_en5141", subset = "Fuel use in different heating types") # [[Energy balance]]
colnames(fuelSharesgeneric@data) <- gsub("[ \\.]", "_", colnames(fuelSharesgeneric@data))
#fuelSharesgeneric@data <- merge(fuelSharesgeneric@data, data.frame(Time = 1900:2080))

objects.store(fuelSharesgeneric)
cat("Object fuelSharesgeneric initiated!\n")

+ Show code - Hide code

## This is code Op_en5141/fuelUse on page [[Energy balance]].
library(OpasnetUtils)

fuelUse <- Ovariable("fuelUse",
	dependencies = data.frame(
		Name = c("energyUse", "fuelShares")
	),
	formula = function(...) {
		out <- energyUse * fuelShares
		return(out)
	}
)

objects.store(fuelUse)
cat("Ovariable fuelUse stored.\n")

Old version with a set of linear equations

Energy balances are described as input = output on a coarse level (called classes) where the structure is the same or similar to the OECD energy balance tables. If possible, this is described on the Energy balance method level and it is shared by all cities.
On more detailed (variable level in the matrix), the fraction of each variable of the total class are described separately. Fractions are city specific and they are described on city level in a separate table.
Based on the fraction table, detailed equations with variables are created. The format will be fraction * class total = variable.
The last fraction has zero degrees of freedom when the class total is given. However, it must have a variable and thus a row in the fraction table. The result for that variable is an empty cell (which results in NA).
Unlike in the previous version, all variables are given either as values or equations, and the user interface is not used for BAU. In contrast, user interface or decision table may be used to derive values for alternative scenarios.
To make this work, the city-specific fraction data must be defined as ovariable (so that it can be changed with a decision table), and also the energy balance method must be described asa ovariable. How are we going to make the two interplay, as we may want to have several cities?
- Define one city ovariable and evaluate energy balance with that. The ovariable has a generic name. Then, define a new city ovariable with the same name and re-evaluate the energy balance ovariable; this must be done so that the two cities are appended rather than replaced.
- city ovariables are appended first into a large fraction table, and then that is used to create the large energy balance matrix. ←--#: . This is clearly better. --Jouni 17:09, 21 February 2013 (EET) (type: truth; paradigms: science: defence)
The city-specific ovariable may have Iter and other indices. A separate matrix is created and solved for each unique combination of indices. This makes it possible to have a very flexible approach.
We should check if the energy balance matrix (see Matti's Excel) has city-specific equations. If possible, energy transformations are described as generic equations on the energy balance method.
Structure of OECD Energy balance tables (data):
- Fuel (given as observation columns in OECD table)
- Activity (row in OECD table)
- Description
Structure of the generic process table
- Equation,
- Col,
- Result,
- Description? ⇤--#: . This does not join up in a coherent way. --Jouni 17:09, 21 February 2013 (EET) (type: truth; paradigms: science: attack)
Columns for fraction table
- Class
- Item
- Result (fraction)
- Indices as needed

Example table for making matrices from text format equations. CHPcapacity describes which of the piecewise linear equations should be used. Policy is a decision option that alters the outcome. Dummy is only for compatibility but it is not used.

Equations(GWh /a)
Obs	CHPcapacity	Policy	Equation	Dummy	Description
1		Biofuel	CHP renewable = CHP peat	1	Biofuel policy contains half biofuels, half peat
2		BAU	CHP renewable = 89.24	1
3			CHP peat + CHP renewable + CHP oil = CHP heat + CHP electricity	1
4			CHP peat = 90-98*CHP oil	1
5			CHP electricity = 0.689*CHP heat	1
6	CHP<1000		H heat = 0.08*CHP heat	1	Small heat plants reflect the total heat need
7	CHP>1000		CHP heat + CHP electricity = 1000	1	But production capacity of CHP may be overwhelmed, decoupling CHP heat and H heat.
8			H biogas + H oil = H heat	1
9			H oil = 18.973*H biogas	1
10			Bought electricity + CHP electricity = Cons electricity	1
11			CHP heat + H heat = Cons heat	1
12			Cons electricity = 900-1100	1
13			Cons heat = 900-1000	1

Example table to describe the details about nonlinear equations.

Nonlinearity parameters(GWh /a)
Obs	critVar	critIndex	rescol	critLocLow	critLocHigh	critValue
1	Cons heat	CHPcapacity	Result	CHP<1000	CHP>1000	1080

This table is fetched if there are no nonlinearities. Therefore, there is no need to copy it to the case study page.

No nonlinearities(GWh /a)
Obs	critVar	critIndex	rescol	critLocLow	critLocHigh	critValue
1

This table is fetched if there are no modelled upstream variables that would affect the equations.

No modelled upstream variables(-)
Obs	energybalanceVars	Result
1

Stored objects below used by Energy balance in Kuopio.

+ Show code - Hide code

library(OpasnetUtils)
library(ggplot2)

solveMatrix <- function(equations, directinputtemp, N. = 0, verbose = FALSE, dbug = FALSE, ...) {
	if (N. == 0) N. <- openv$N
	N <- N.
	
	# Parse equations using regular expressions
	
	out <- gsub(" ", "", equations)
	out <- strsplit(out, "=")
	out <- lapply(out, strsplit, split = "+", fixed = TRUE)
	nterms <- rapply(out, length) # number of terms on each side of each equation
	signs <- rep(
			rep(c(1,-1), length(out)), # lhs - rhs = 0
			nterms
	)
	eqn <- factor(rep(rep(1:length(equations), each = 2), nterms)) # equation number
	out <- unlist(out)
	k <- rep(1, length(out))
	x <- out
	
	# Equation factors
	
	out <- strsplit(out, split = "*", fixed = TRUE)
	temp <- sapply(out, length) > 1
	
	for (i in (1:length(out))[temp]) {
		k[i] <- out[[i]][1]
		x[i] <- out[[i]][2]
	}
	
	# Constants
	
	temp <- substr(x, 1, 1) %in% 0:9 # is interpretable (probably)
	
	for (i in (1:length(out))[temp]) {
		k[i] <- x[i]
		x[i] <- "Constant"
	}
	
	# Direct input (substitutes for constants)
	
	if (nrow(directinputtemp) > 0) {
		l <- tapply(eqn, data.frame(eqn), length)
		di <- directinputtemp$energybalanceVars
		dires <- directinputtemp$directinputResult
		for (i in 1:length(di)) {
			temp <- eqn %in% names(l[l == 2]) & x %in% c("Constant", di[i])
			temp <- tapply(temp, eqn, sum)
			k[eqn %in% names(temp[temp == 2]) & x == "Constant"] <- dires[i]
		}
	}
	
	# Interpret possible distribution regular expressions (also converts text to numeric)
	
	out <- data.frame(eqn, x, signs, Result = k)
	out <- interpret(out, N = N)
	
	marg <- colnames(out)[colnames(out) %in% c("eqn", "x", "Iter")]
	
	out$k <- out$Result * out$signs
	
	# Separate constants and vars and make them into arrays
	
	temp <- out[out$x != "Constant", ]
	
	b <- out[out$x == "Constant", ]
	b <- tapply(b$k, b[marg[-2]], I)
	b[is.na(b)] <- 0
	
	temp <- tapply(temp$k, temp[marg], I)
	if ("Iter" %in% marg) {
		temp <- temp[,dimnames(temp)$x != "Constant",]
	} else {
		temp <- temp[,dimnames(temp)$x != "Constant"]
	}
	temp[is.na(temp)] <- 0
	
	if ("Iter" %in% marg) {
		out <- list()
		for (i in 1:N) {
			ret <- tryCatch(
				out[[i]] <- solve(temp[,,i], -b[,i]),
				error = function(...) return(NULL)
			)
			if (is.null(ret)) {
				print("Faulty matrix:")
				oprint(temp[,,i])
				stop(geterrmessage())
			}
		}
		out <- data.frame(
			energybalanceVars = names(unlist(out)), 
			Iter = rep(1:N, each = length(out[[1]])), 
			Result = unlist(out)
		)
	} else {
		ret <- tryCatch(
			out <- solve(temp, -b),
			error = function(...) return(NULL)
		)
		if (is.null(ret)) {
			print("Faulty matrix:")
			oprint(temp)
			stop(geterrmessage())
		}
		out <- data.frame(
			energybalanceVars = names(out),
			Result = out
		)
	}
	return(out)
}

energy.balance <- Ovariable( # NOTE! energy.balance requires that there exists object N that defines the number of iterations.
	name = "energy.balance",
	dependencies = data.frame(
		Name = c("balance", "nonlinearity", "directinput")
	),
	formula = function(...) {
		# Column Equation: change factor into character and demarginalise.
		balance@data$Equation <- as.character(balance@data$Equation)
		balance@marginal[colnames(balance@output) == "Equation"] <- FALSE
		
		# Take those indices that are not needed within solveMatrix and merge them to the balance from the directinput.
		directinputindex <- directinput@output[directinput@marginal]
		directinputindex <- directinputindex[!colnames(directinputindex) %in% c("Iter", "energybalanceVars")]
		
		indices <- unique(merge(balance@output[balance@marginal], directinputindex))
		
		# Look at each unique combination of index locations and solve the set of equations to each.
		
		out <- data.frame()
		for(k in 1:nrow(indices)) {
			
			equations <- merge(balance@output, indices[k , , drop = FALSE])
			
			directinputtemp <- merge(directinput@output[!is.na(result(directinput)),], indices[k, colnames(indices) != "Iter", drop = FALSE])
			directinputtemp <- directinputtemp[ , c("energybalanceVars", "directinputResult")]
			
			temp <- solveMatrix(as.character(equations[ , "Equation"]), directinputtemp, N. = N, ...)
			
			out <- rbind(out, merge(indices[k , , drop = FALSE], temp))
		}
		
		if(!(nonlinearity@output$critVar[1] == "" | is.na(nonlinearity@output$critVar[1]))) {
			
			# Go through the non-linear equations one at a time and find which values to use in which situation.
			
			for(i in 1:nrow(nonlinearity@output)) {
				
				critVar 	<- gsub(" ", "", as.character(nonlinearity@output$critVar[i]))
				critIndex 	<- as.character(nonlinearity@output$critIndex[i])
				rescol 		<- as.character(nonlinearity@output$rescol[i])
				critLocLow 	<- as.character(nonlinearity@output$critLocLow[i])
				critLocHigh	<- as.character(nonlinearity@output$critLocHigh[i])
				critValue 	<- result(nonlinearity)[i]
				
				marginal <- colnames(out) %in% c(colnames(balance@output)[balance@marginal], rescol, "Iter")
				
				choos <- out[out$energybalanceVars == critVar , marginal]
				choos <- choos[
					choos[critIndex] == critLocLow  & choos[rescol] <= critValue | 
						choos[critIndex] == critLocHigh & choos[rescol] >  critValue , 
					colnames(choos) != rescol # Remove result column 
				]
				
				out <- merge(out, choos)
			}
		}
		
		# Create the marginal: those from the parents and "energybalanceVars" but not nonlinear indices.
		outmarginal <- c(
			colnames(balance@output)[balance@marginal], 
			colnames(nonlinearity@output)[nonlinearity@marginal], 
			colnames(directinput@output)[directinput@marginal],
			"energybalanceVars"
		)
		outmarginal <- outmarginal[!outmarginal %in% nonlinearity@output$critIndex]
		outmarginal <- colnames(out) %in% outmarginal
		
		out <- new("ovariable", name = "energy.balance", output = out, marginal = outmarginal)
		
		return(out)
	}
)

objects.store(solveMatrix, energy.balance)

cat("Objects solveMatrix and energy.balance stored with code initiate.\n")

Model version that was used to run results for ISEE2013.

How to give uncertain parameters?

In equations, the content is interpreted only inside solveMatrix. Therefore, the typical approach where all unique index combinations are run one at a time does not work.
There should be an update in parameter interpretation for terms with one entry only. It can no longer be based on as.numeric, if distributions (=text) is allowed.
- If it starts with [a-z.] it is a variable name.
- If it starts with [0-9<\\-] it is a parameter value.
Instead of params[[i]] and [[vars]] vectors, a data.frame will be created with Result as the params column.
The data.frame is then interpreted with N = N. If parameters are probabilistic, Iter column will appear.
When all parameters have been interpreted, check if Iter exists.
If Iter exists, make a for loop for all values of Iter.
- Create a matrix from the parameters and solve.
- Rbind the result to a data.frame with Iter.
Return the output.

Old code with an input table with columns Equation, Col, Result, Description: [1]

References

Related files

In English
Assessment	Main page \| Helsinki energy decision options 2015
Helsinki data	Building stock in Helsinki \| Helsinki energy production \| Helsinki energy consumption \| Energy use of buildings \| Emission factors for burning processes \| Prices of fuels in heat production \| External cost
Models	Building model \| Energy balance \| Health impact assessment \| Economic impacts
Related assessments	Climate change policies in Helsinki \| Climate change policies and health in Kuopio \| Climate change policies in Basel
In Finnish
Yhteenveto	Helsingin energiapäätös 2015 \| Helsingin energiapäätöksen vaihtoehdot 2015 \| Helsingin energiapäätökseen liittyviä arvoja \| Helsingin energiapäätös 2015.pptx

Energy balance: Difference between revisions

Latest revision as of 07:49, 10 April 2019

Contents

Question

Answer

Rationale

Energy balance with linear programming

Fuel use and fuel shares in generic processes

Old version with a set of linear equations

See also

References

Related files

Navigation menu

@@ Line 1: / Line 1: @@
+<noinclude>
 [[op_fi:Energiatase]]
 [[Category:Energy]]
 [[Category:Energy balance]]
 [[Category:Online model]]
+[[Category:Contains R code]]
+[[Category:Code under inspection]]
 {{method|moderator=Jouni|stub=Yes}}
+</noinclude>
-{{attack|# |Don't put arguments on the same line as a subheading. Check the subheadings, you have some that are not used in [[method]]s.|--[[User:Jouni|Jouni]] 22:46, 12 February 2013 (EET)}}
+==Question==
+What is energy balance and how is it modelled?
-The energy calculation based on two core concepts
+==Answer==
-.Energy production or source
+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.
-.Energy consumption
+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]].
+<gallery heights="500" widths="500">
+File:District heat flow in Helsinki.png|Example of energy balance model: District heat flow in Helsinki. The scenarios are from the assessment [[Helsinki energy decision 2015]].
+File:Flow of electric power in Helsinki.png|Example of energy balance model: Electric power heat flow in Helsinki. The scenarios are from the assessment [[Helsinki energy decision 2015]].
+File:Incomes and costs of energy production in Helsinki.png|Example of energy balance model: Incomes and costs of energy production in Helsinki. The scenarios are from the assessment [[Helsinki energy decision 2015]].
+</gallery>
-==Question==
+This code is an example how the energy balance model is used in a city case. The data comes from [[Helsinki energy decision 2015]].
-How to calculate energy balances?
+<rcode name="answer" graphics="1">
+library(OpasnetUtils)
+library(ggplot2)
-{{comment|# |Clear question|--[[User:Sam0911|Sam0911]] 10:28, 17 February 2013 (EET)}}
+openv.setN(0) # use medians instead of whole sampled distributions
-===Purpose===
+# Download case-specific data, in this case from Helsinki.
-where is the purpose?
+objects.latest("Op_en7237", code_name = "intermediates") # [[Helsinki energy decision 2015]]
-===Boundaries===
+# Download energy balance model and its parts:
-no info about bouderies?
+#  EnergyConsumerDemandTotal
+#  EnergyFlowHeat
+#  EnergyFlowOther
+#  EnergyNetworkDemand
+#  EnergyNetworkOptim
+#  fuelUse
+#  EnergyNetworkCost
+objects.latest("Op_en5141", code_name = "EnergyNetworkOptim") # [[Energy balance]]
-===Procedure===
+oprint(summary(EvalOutput(EnergyNetworkOptim)))
-{{attack|# |Where is the procedure of the assessment?|--[[User:Sam0911|Sam0911]] 10:13, 17 February 2013 (EET)}}
-===Management===
+</rcode>
-{{attack|# |Where is the management of the assessment?|--[[User:Sam0911|Sam0911]] 10:13, 17 February 2013 (EET)}}
 ==Rationale==
-===Input===
+=== Energy balance with linear programming ===
-===Output===
+The linear programming problem is set up as follows.
-===Data===
+For each production unit: let x<sub>i</sub> be activity of the plant. Lets also have variables y<sub>j</sub> for deficits and excesses for each type of energy produced.
-===Calculations===
+The objective function is the function we are optimising. Each production unit has a unit profit per activity denoted by a<sub>i</sub> 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 b<sub>j</sub>. The whole objective function then becomes: sum(x<sub>i</sub>a<sub>i</sub>) + sum(y<sub>j</sub>b<sub>j</sub>).
-===Rationale===
-{{attack|# |What are other energy sources,i wish it should be specific,how about air transport,does it included in the energy balance calculation?|--[[User:Sam0911|Sam0911]] 16:29, 10 February 2013 (EET)}}
+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.
+:''Ovariables like EnergyNetworkOptim below are used in [[Helsinki energy decision 2015]]. [[Prices of fuels in heat production]] are used as direct inputs in the optimising.
+<rcode name="EnergyNetworkOptim" label="Initiate EnergyNetworkOptim(developers only)" embed=1>
+## This code is  Op_en5141/EnergyNetworkOptim [[Energy balance]].
+library(OpasnetUtils)
-Use a table where different fuel types are columns and different stocks, energy production processes, or consumption types are rows called ''Energy accounts'' (Ene.account for short). See also an example [[File:Energy supply in Europe.xls]]. All Ene.accounts and fuel types used are listed on [[Energy consumption classes]].
+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)
+	}
+)
-{| {{prettytable}}
+EnergyFlowHeat <- Ovariable("EnergyFlowHeat",
-|+ '''A part of a full energy balance sheet (ktoe/a).
+	dependencies = data.frame(
-|----
+		Name = c(
-!
+			"EnergyConsumerDemandTotal",
-! Coal and peat
+			"fuelShares"
-! Crude oil
+		)
-! Petrochemical products
+	),
-! Electricity
+	formula = function(...) {
-! Heat
+		EnergyFlowHeat <- EnergyConsumerDemandTotal[EnergyConsumerDemandTotal$Consumable %in% c("Heating", "Hot water"),] * fuelShares
-|----
+		return(EnergyFlowHeat)
-!colspan="6"|Production and supply of primary energy
+	}
-|----
+)
-| Production || 129|| || || ||
-|----
-| Import || || 28|| 63|| ||
-|----
-| Export || || || || ||
-|----
-!colspan="6"| Conversion of primary energy to use
-|----
-| Transfers || || || || ||
-|----
-| Electricity plants || || || || ||
-|----
-| CHP plants || -129|| -28|| || 61|| 96
-|----
-!colspan="6"|Final energy consumption
-|----
-| Industry || || || || ||
-|----
-| Road transport || || || 54|| ||
-|----
-| Heating || || || || 10|| 96
-|----
-| Other energy use || || || 9|| 51||
-|}
-Energy balance can be considered as double-entry bookkeeping. The production and conversion are typically credit (i.e., where the energy comes from), while consumption is debit (i.e., where the energy is used). In a typical case, both credit and debit are marked positive on the energy balance sheet. However, there are activities that convert one fuel to another, so that the energy is moved from one column to another. In this case, credit (the source) is marked negative and debit (the target) is marked positive. This is because the production and conversion rows together should reflect how much energy is actually available for final consumption. In other words, the sum of production and conversion rows should equal the sum of the consumption rows in every column.
+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)
+	}
+)
-Because production + conversion = consumption, also production = consumption - conversion. These equations are used to derive the supply that is needed to fulfil the demand.
+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)
+	}
+)
-{{comment|# |it seems the formula is ok,the code also runs|--[[User:Sam0911|Sam0911]] 16:33, 10 February 2013 (EET)}}=== Formula ===
+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)
-Default run key: qffJP3u5RpFw1GuL
+		#out <- data.frame()
-<rcode variables="
+		include <- union(ePmargs, pPmargs)
-	name:Frac_CHP_renew|description:Amount of renewable fuels compared with peat (may be over 1)|default:0.0544|
+		include <- union(include, ENDmargs)
-	name:Ind_process_heat|description:Consumption of process heat produced by industry itself|default:685.4|
+		include <- union(include, fPmargs)
-	name:Cons_Home_electricity|description:Consumption of electricity by households|default:296|
+		include <- include[!include %in% exclude]
-	name:Cons_Ind_electricity|description:Consumption of electricity by industry|default:204|
-	name:Cons_Commerce_electricity|description:Consumption of electricity by commerce|default:334|
+		# Optimisation iteration function
-	name:Cons_Municip_electricity|description:Consumption of electricity by municipality|default:50.3|
-	name:Cons_Home_heat|description:Consumption of district heat by households|default:303|
+		optf <- function(varsi) {
-	name:Cons_Ind_heat|description:Consumption of district heat by industry|default:57.1|
-	name:Cons_Commerce_heat|description:Consumption of district heat by commerce|default:224|
+			# Acquire relevant sections of variables with respect to index iteration
-	name:Cons_Municip_heat|description:Consumption of district heat by municipality|default:386
-">
+			demandi <- unlist(varsi[1,grepl("Demand$", colnames(varsi))])
-library(OpasnetUtils)
+			names(demandi) <- gsub("Demand$", "", names(demandi))
-library(xtable)
+			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(
+,
+				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)
-# These three are city-specific: Matrix_page, Additions, Solutions.
+				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)),
+:length(plants)
+			]) <- 1
+			diag(M[
+				(rowi + length(plants) + 1):(rowi + 2 * length(plants)),
+: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) +
+:(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
-Matrix_page <- "Op_en5469" # This is the page with the matrix content.
+			# 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)))
-Additions <- data.frame( # These are additional contents to the matrix, in addition to Matrix_page.
+		return(out)
-	Equation = "Fraction CHP renewable",
+	}
-	Col = "CHP peat",
-	ETMResult = Frac_CHP_renew,
-	ETMDescription = "Actual number given by user."
 )
-Solutions <- c(rep(0, 18), # These are the solutions (right sides) of the equations.
+fuelUse <- Ovariable("fuelUse",
-	Ind_process_heat,
+	dependencies = data.frame(
-	Cons_Home_electricity,
+		Name = c(
-	Cons_Ind_electricity,
+			"EnergyNetworkOptim",
-	Cons_Commerce_electricity,
+			"EnergyNetworkDemand",
-	Cons_Municip_electricity,
+			"energyProcess",
-	Cons_Home_heat,
+			"EnergyFlowHeat",
-	Cons_Ind_heat,
+			"temperdays"
-	Cons_Commerce_heat,
+		)
-	Cons_Municip_heat
+	),
+	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)
+	}
 )
-dependencies <- data.frame(Name = c(
+EnergyNetworkCost <- Ovariable("EnergyNetworkCost",
-	"Matrix_page",
+	dependencies =
-	"Additions",
+		data.frame(
-	"Solutions"
+			Name = c(
-), Key = "")
+				"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)
-formula <- function(dependencies, intermediates = FALSE, ...) {
+		maintenance <- oper[oper@output$Parameter == "Management cost" , colnames(oper@output) != "Parameter"]
+		maintenance <- unkeep(maintenance, sources = TRUE)
+		#maintenance <- oapply(maintenance, cols = c("Plant"), FUN = sum)
-	ComputeDependencies(dependencies, ...)
+		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)
+	}
+)
-	## Define the energy transformation matrix ETM.
+objects.store(EnergyConsumerDemandTotal, EnergyFlowHeat, EnergyFlowOther, EnergyNetworkDemand, EnergyNetworkOptim, fuelUse, EnergyNetworkCost)
+cat("Ovariables EnergyConsumerDemandTotal, EnergyFlowHeat, EnergyFlowOther, EnergyNetworkDemand, EnergyNetworkOptim, EnergyNetworkCost and fuelUse stored.\n")
+</rcode>
-	# Get the data from Opasnet Base.
+=== Fuel use and fuel shares in generic processes ===
-	ETM <- tidy(op_baseGetData("opasnet_base", Matrix_page), objname = "ETM", direction = "wide")
+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.
-	# Add user-defined data to it.
+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.
-	ETM <- rbind(ETM, Additions)
+<t2b name='Fuel use in different heating types' index='Heating,Burner,Fuel' obs='Fraction' desc='Description' unit='-'>
+Wood|Household|Wood|1|
+Oil|Household|Light oil|1|
+Gas|Household|Gas|1|
+Heating oil|Household|Light oil|1|
+Coal|Household|Coal|1|
+Other sources|Household|Other sources|1|
+No energy source|Household|Other sources|1|
+Geothermal|Grid|Electricity|0.3|Geothermal does not sum up to 1 because more heat is produced than electricity consumed.
+Centrifuge, hydro-extractor|Grid|Electricity|0.3|Not quite clear what this is but presumably a heat pump.
+Solar heater/ collector|Grid|Electricity|0.1|Use only; life-cycle impacts omitted.
+Electricity|Grid|Electricity|1|
+District|Undefined|Heat|1|
+</t2b>
-	# Reshape the table to make a matrix, and replace NA with 0.
+<rcode name='fuelSharesgeneric' embed=1 label='Initiate fuelSharesgeneric (only for developers)'>
-	ETM <- reshape(ETM, timevar = "Col", idvar = "Equation", v.names = "ETMResult", direction = "wide")
+# This is code Op_en5141/fuelSharesgeneric (only generic) on page [[Energy balance]].
-	colnames(ETM) <- gsub("ETMResult.", "", colnames(ETM))
+library(OpasnetUtils)
-	if(intermediates) {print(xtable(ETM), type = 'html')} # Show the matrix with explanations.
+fuelSharesgeneric <- Ovariable("fuelSharesgeneric", ddata = "Op_en5141", subset = "Fuel use in different heating types") # [[Energy balance]]
+colnames(fuelSharesgeneric@data) <- gsub("[ \\.]", "_", colnames(fuelSharesgeneric@data))
+#fuelSharesgeneric@data <- merge(fuelSharesgeneric@data, data.frame(Time = 1900:2080))
-	ETM <- as.matrix(ETM[, !colnames(ETM) %in% c("Equation", "ETMDescription")])
+objects.store(fuelSharesgeneric)
+cat("Object fuelSharesgeneric initiated!\n")
-	for(i in 1:ncol(ETM)) {
+</rcode>
-		ETM[is.na(ETM[, i]), i] <- 0
+<rcode name="fuelUse" label="Initiate fuelUse (developers only)" embed=1>
+## This is code Op_en5141/fuelUse on page [[Energy balance]].
+library(OpasnetUtils)
+fuelUse <- Ovariable("fuelUse",
+	dependencies = data.frame(
+		Name = c("energyUse", "fuelShares")
+	),
+	formula = function(...) {
+		out <- energyUse * fuelShares
+		return(out)
 	}
+)
+objects.store(fuelUse)
+cat("Ovariable fuelUse stored.\n")
+</rcode>
+<noinclude>
+=== Old version with a set of linear equations ===
+* Energy balances are described as input = output on a coarse level (called classes) where the structure is the same or similar to the OECD energy balance tables. If possible, this is described on the Energy balance method level and it is shared by all cities.
+* On more detailed (variable level in the matrix), the fraction of each variable of the total class are described separately. Fractions are city specific and they are described on city level in a separate table.
+* Based on the fraction table, detailed equations with variables are created. The format will be fraction * class total = variable.
+* The last fraction has zero degrees of freedom when the class total is given. However, it must have a variable and thus a row in the fraction table. The result for that variable is an empty cell (which results in NA).
+* Unlike in the previous version, all variables are given either as values or equations, and the user interface is not used for BAU. In contrast, user interface or decision table may be used to derive values for alternative scenarios.
+* To make this work, the city-specific fraction data must be defined as ovariable (so that it can be changed with a decision table), and also the energy balance method must be described asa ovariable. How are we going to make the two interplay, as we may want to have several cities?
+** Define one city ovariable and evaluate energy balance with that. The ovariable has a generic name. Then, define a new city ovariable with the same name and re-evaluate the energy balance ovariable; this must be done so that the two cities are appended rather than replaced.
+** city ovariables are appended first into a large fraction table, and then that is used to create the large energy balance matrix. {{defend|# |This is clearly better.|--[[User:Jouni|Jouni]] 17:09, 21 February 2013 (EET)}}
+* The city-specific ovariable may have Iter and other indices. A separate matrix is created and solved for each unique combination of indices. This makes it possible to have a very flexible approach.
+* We should check if the energy balance matrix (see Matti's Excel) has city-specific equations. If possible, energy transformations are described as generic equations on the energy balance method.
+* Structure of OECD Energy balance tables (data):
+** Fuel (given as observation columns in OECD table)
+** Activity (row in OECD table)
+** Description
+* Structure of the generic process table
+** Equation,
+** Col,
+** Result,
+** Description? {{attack|# |This does not join up in a coherent way.|--[[User:Jouni|Jouni]] 17:09, 21 February 2013 (EET)}}
+* Columns for fraction table
+** Class
+** Item
+** Result (fraction)
+** Indices as needed
+; Example table for making matrices from text format equations. CHPcapacity describes which of the piecewise linear equations should be used. Policy is a decision option that alters the outcome. Dummy is only for compatibility but it is not used.
+<t2b name="Equations" index="CHPcapacity,Policy,Equation" obs="Dummy" desc="Description" unit="GWh /a">
+|Biofuel|CHP renewable = CHP peat|1|Biofuel policy contains half biofuels, half peat
+|BAU|CHP renewable = 89.24|1|
+||CHP peat + CHP renewable + CHP oil = CHP heat + CHP electricity|1|
+||CHP peat = 90-98*CHP oil|1|
+||CHP electricity = 0.689*CHP heat|1|
+CHP<1000||H heat = 0.08*CHP heat|1|Small heat plants reflect the total heat need
+CHP>1000||CHP heat + CHP electricity = 1000|1|But production capacity of CHP may be overwhelmed, decoupling CHP heat and H heat.
+||H biogas + H oil = H heat|1|
+||H oil = 18.973*H biogas|1|
+||Bought electricity + CHP electricity = Cons electricity|1|
+||CHP heat + H heat = Cons heat|1|
+||Cons electricity = 900-1100|1|
+||Cons heat = 900-1000|1|
+</t2b>
+; Example table to describe the details about nonlinear equations.
+<t2b name="Nonlinearity parameters" index="critVar,critIndex,rescol,critLocLow,critLocHigh" obs="critValue" unit="GWh /a">
+Cons heat|CHPcapacity|Result|CHP<1000|CHP>1000|1080
+</t2b>
+; This table is fetched if there are no nonlinearities. Therefore, there is no need to copy it to the case study page.
-	# Define the outputs for equations (for most equations, it is a zero sum game).
+<t2b name="No nonlinearities" index="critVar,critIndex,rescol,critLocLow,critLocHigh" obs="critValue" unit="GWh /a">
+|||||
+</t2b>
-	out <- solve(ETM, Solutions) # Solve the matrix.
+; This table is fetched if there are no modelled upstream variables that would affect the equations.
-	out <- data.frame(Energy.class = colnames(ETM), Result = out)
+<t2b name="No modelled upstream variables" index="energybalanceVars" obs="Result" unit="-">
+|
+</t2b>
+Stored objects below used by [[Energy balance in Kuopio]].
+<rcode name="initiate" label="Initiate objects" store=1>
+library(OpasnetUtils)
+library(ggplot2)
+solveMatrix <- function(equations, directinputtemp, N. = 0, verbose = FALSE, dbug = FALSE, ...) {
+	if (N. == 0) N. <- openv$N
+	N <- N.
+	# Parse equations using regular expressions
+	out <- gsub(" ", "", equations)
+	out <- strsplit(out, "=")
+	out <- lapply(out, strsplit, split = "+", fixed = TRUE)
+	nterms <- rapply(out, length) # number of terms on each side of each equation
+	signs <- rep(
+			rep(c(1,-1), length(out)), # lhs - rhs = 0
+			nterms
+	)
+	eqn <- factor(rep(rep(1:length(equations), each = 2), nterms)) # equation number
+	out <- unlist(out)
+	k <- rep(1, length(out))
+	x <- out
+	# Equation factors
+	out <- strsplit(out, split = "*", fixed = TRUE)
+	temp <- sapply(out, length) > 1
+	for (i in (1:length(out))[temp]) {
+		k[i] <- out[[i]][1]
+		x[i] <- out[[i]][2]
+	}
+	# Constants
+	temp <- substr(x, 1, 1) %in% 0:9 # is interpretable (probably)
+	for (i in (1:length(out))[temp]) {
+		k[i] <- x[i]
+		x[i] <- "Constant"
+	}
+	# Direct input (substitutes for constants)
+	if (nrow(directinputtemp) > 0) {
+		l <- tapply(eqn, data.frame(eqn), length)
+		di <- directinputtemp$energybalanceVars
+		dires <- directinputtemp$directinputResult
+		for (i in 1:length(di)) {
+			temp <- eqn %in% names(l[l == 2]) & x %in% c("Constant", di[i])
+			temp <- tapply(temp, eqn, sum)
+			k[eqn %in% names(temp[temp == 2]) & x == "Constant"] <- dires[i]
+		}
+	}
+	# Interpret possible distribution regular expressions (also converts text to numeric)
+	out <- data.frame(eqn, x, signs, Result = k)
+	out <- interpret(out, N = N)
+	marg <- colnames(out)[colnames(out) %in% c("eqn", "x", "Iter")]
+	out$k <- out$Result * out$signs
+	# Separate constants and vars and make them into arrays
+	temp <- out[out$x != "Constant", ]
+	b <- out[out$x == "Constant", ]
+	b <- tapply(b$k, b[marg[-2]], I)
+	b[is.na(b)] <- 0
+	temp <- tapply(temp$k, temp[marg], I)
+	if ("Iter" %in% marg) {
+		temp <- temp[,dimnames(temp)$x != "Constant",]
+	} else {
+		temp <- temp[,dimnames(temp)$x != "Constant"]
+	}
+	temp[is.na(temp)] <- 0
+	if ("Iter" %in% marg) {
+		out <- list()
+		for (i in 1:N) {
+			ret <- tryCatch(
+				out[[i]] <- solve(temp[,,i], -b[,i]),
+				error = function(...) return(NULL)
+			)
+			if (is.null(ret)) {
+				print("Faulty matrix:")
+				oprint(temp[,,i])
+				stop(geterrmessage())
+			}
+		}
+		out <- data.frame(
+			energybalanceVars = names(unlist(out)),
+			Iter = rep(1:N, each = length(out[[1]])),
+			Result = unlist(out)
+		)
+	} else {
+		ret <- tryCatch(
+			out <- solve(temp, -b),
+			error = function(...) return(NULL)
+		)
+		if (is.null(ret)) {
+			print("Faulty matrix:")
+			oprint(temp)
+			stop(geterrmessage())
+		}
+		out <- data.frame(
+			energybalanceVars = names(out),
+			Result = out
+		)
+	}
 	return(out)
 }
-energy.balance <- new("ovariable",
+energy.balance <- Ovariable( # NOTE! energy.balance requires that there exists object N that defines the number of iterations.
 	name = "energy.balance",
-	dependencies = dependencies,
+	dependencies = data.frame(
-	formula = formula
+		Name = c("balance", "nonlinearity", "directinput")
+	),
+	formula = function(...) {
+		# Column Equation: change factor into character and demarginalise.
+		balance@data$Equation <- as.character(balance@data$Equation)
+		balance@marginal[colnames(balance@output) == "Equation"] <- FALSE
+		# Take those indices that are not needed within solveMatrix and merge them to the balance from the directinput.
+		directinputindex <- directinput@output[directinput@marginal]
+		directinputindex <- directinputindex[!colnames(directinputindex) %in% c("Iter", "energybalanceVars")]
+		indices <- unique(merge(balance@output[balance@marginal], directinputindex))
+		# Look at each unique combination of index locations and solve the set of equations to each.
+		out <- data.frame()
+		for(k in 1:nrow(indices)) {
+			equations <- merge(balance@output, indices[k , , drop = FALSE])
+			directinputtemp <- merge(directinput@output[!is.na(result(directinput)),], indices[k, colnames(indices) != "Iter", drop = FALSE])
+			directinputtemp <- directinputtemp[ , c("energybalanceVars", "directinputResult")]
+			temp <- solveMatrix(as.character(equations[ , "Equation"]), directinputtemp, N. = N, ...)
+			out <- rbind(out, merge(indices[k , , drop = FALSE], temp))
+		}
+		if(!(nonlinearity@output$critVar[1] == "" | is.na(nonlinearity@output$critVar[1]))) {
+			# Go through the non-linear equations one at a time and find which values to use in which situation.
+			for(i in 1:nrow(nonlinearity@output)) {
+				critVar 	<- gsub(" ", "", as.character(nonlinearity@output$critVar[i]))
+				critIndex 	<- as.character(nonlinearity@output$critIndex[i])
+				rescol 		<- as.character(nonlinearity@output$rescol[i])
+				critLocLow 	<- as.character(nonlinearity@output$critLocLow[i])
+				critLocHigh	<- as.character(nonlinearity@output$critLocHigh[i])
+				critValue 	<- result(nonlinearity)[i]
+				marginal <- colnames(out) %in% c(colnames(balance@output)[balance@marginal], rescol, "Iter")
+				choos <- out[out$energybalanceVars == critVar , marginal]
+				choos <- choos[
+					choos[critIndex] == critLocLow  & choos[rescol] <= critValue |
+						choos[critIndex] == critLocHigh & choos[rescol] >  critValue ,
+					colnames(choos) != rescol # Remove result column
+				]
+				out <- merge(out, choos)
+			}
+		}
+		# Create the marginal: those from the parents and "energybalanceVars" but not nonlinear indices.
+		outmarginal <- c(
+			colnames(balance@output)[balance@marginal],
+			colnames(nonlinearity@output)[nonlinearity@marginal],
+			colnames(directinput@output)[directinput@marginal],
+			"energybalanceVars"
+		)
+		outmarginal <- outmarginal[!outmarginal %in% nonlinearity@output$critIndex]
+		outmarginal <- colnames(out) %in% outmarginal
+		out <- new("ovariable", name = "energy.balance", output = out, marginal = outmarginal)
+		return(out)
+	}
 )
-temp <- EvalOutput(energy.balance)
+objects.store(solveMatrix, energy.balance)
-print(xtable(temp@output), type = 'html')
-objects.put(energy.balance)
+cat("Objects solveMatrix and energy.balance stored with code initiate.\n")
-cat("Object energy.balance saved. Write down the key of the model run page for further use.\n")
+</rcode>
-</rcode>
-The model first creates the matrix (and shows it to the user if the user has defined intermediates = TRUE), and then it fills all empty cells with 0 and solves it by using the input values the user has given via the user interface. (The user interface can later be replaced by another model.) All the current input values are final consumptions. The example interface and values are based on the Kuopio data, but they are not stored as a part of the ovariable.
+* [http://en.opasnet.org/en-opwiki/index.php?title=Energy_balance&oldid=30561 Model version] that was used to run results for ISEE2013.
-{{comment|# |good to have a reference to compare|--[[User:Sam0911|Sam0911]] 16:12, 10 February 2013 (EET)}}==See also==
+How to give uncertain parameters?
+* In equations, the content is interpreted only inside solveMatrix. Therefore, the typical approach where all unique index combinations are run one at a time does not work.
+* There should be an update in parameter interpretation for terms with one entry only. It can no longer be based on as.numeric, if distributions (=text) is allowed.
+** If it starts with [a-z.] it is a variable name.
+** If it starts with [0-9<\\-] it is a parameter value.
+* Instead of <nowiki>params[[i]] and [[vars]]</nowiki> vectors, a data.frame will be created with Result as the params column.
+* The data.frame is then interpreted with N = N. If parameters are probabilistic, Iter column will appear.
+* When all parameters have been interpreted, check if Iter exists.
+* If Iter exists, make a for loop for all values of Iter.
+** Create a matrix from the parameters and solve.
+** Rbind the result to a data.frame with Iter.
+* Return the output.
+* Old code with an input table with columns Equation, Col, Result, Description: [http://en.opasnet.org/en-opwiki/index.php?title=Energy_balance&oldid=30515#Code]
+==See also==
+{{Helsinki energy decision 2015}}
+* A previous model version written by Jouni. It used linear optimising but was not fully developed and was not actually used in any final assessment. [http://en.opasnet.org/en-opwiki/index.php?title=Energy_balance&oldid=38536#Calculations]
+* A new version of energy balance is using [https://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html optim function]. See [[Energy balance in Helsinki]].
+* [http://www.nzz.ch/meinung/kommentare/energiepolitik-bitte-wenden-1.18451262 Energiepolitik, bitte wenden!] Neue Zürcher Zeitung 27.12.2014 by Giorgio V. Müller
+* [http://time.com/2981460/electric-cars-home-solar/ TIME: Electric Cars Will Change the Way You Power Your Home]. How the homes of the future will generate and store their own electricity, turning your house into a mini-power plant.
+** [http://time.com/2865496/tesla-opens-patents/ TIME: Tesla opens patents to everyone]
+* [[Energy supply in Europe]]
 * [[media:Health impacts of energy production.ppt]] (a lecture that also contains explanation of an energy balance using matrices)
 * [[Energy balance in Kuopio]] Describes the production and consumption of energy in Kuopio.
+* [[Energy balance in Basel]]
 * [[Energy balance in Stuttgart]] Describes the production and consumption of energy in Stuttgart.
 * [[Energy balance in Suzhou]] Describes the production and consumption of energy in Suzhou.
@@ Line 206: / Line 896: @@
 * [http://www.ymparisto.fi/default.asp?contentid=408832&lan=fi&clan=fi Uusiutuvan energian riskit selvitetään]
 * [http://www.dconnolly.net/research/planning/tools/mesap_planet.html Urgenche: Mesap Planet energy model]
+* [[OpasnetUtils/Drafts]]
 {{urgenche}}
+==References==
+<references/>
 ==Related files==
+</noinclude>
-{{mfiles}}