Difference between revisions of "Energy balance"

library(OpasnetUtils)
library(ggplot2)

objects.latest("Op_en5141", code_name = "initiate")

N <- 10

balance <- Ovariable("balance", ddata = "Op_en5141.equations")

balance@data$Equation <- as.character(levels(balance@data$Equation)[balance@data$Equation]) # Change factor into character.
balance@data$Policy[balance@data$Policy == ""] <- NA # Prepare indices for fillna. This should be part of fillna.
balance@data$CHPcapacity[balance@data$CHPcapacity == ""] <- NA
balance@data <- fillna(balance@data, marginals = 1:2) # Fill empty slots in indices.

nonlinearity <- Ovariable("nonlinearity", ddata = "Op_en5141", subset = "Nonlinearity parameters")

directinput <- Ovariable("directinput", ddata = "Op_en5141", subset = "No modelled upstream variables")


eb <- EvalOutput(energy.balance) 

eb@marginal[colnames(eb@output) == "energy.balanceVar"] <- TRUE # This has to be done manually because CheckMarginals does not notice this as a marginal.
eb@marginal[colnames(eb@output) %in% nonlinearity@output$critIndex] <- FALSE # Nonlinear indices are demarginalised because only one of the two equations apply.

oprint(summary(eb))

ggplot(eb@output, aes_string(x = "energy.balanceVar", y = "energy.balanceResult", fill = "Policy")) + 
	geom_boxplot() +
	theme_grey(base_size = 24) + 
	theme(axis.text.x = element_text(angle = 90, hjust = 1))

# 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")

## This code is  Op_en5141/energyBalance [[Energy balance]].
library(OpasnetUtils)

energyBalance <- Ovariable("energyBalance",
	dependencies = data.frame(
		Name = c(
			"energy", # Energy supply and demand in the system
				#demand and supply that is not centrally controlled.
			"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 outside the balanced processes
			"energyProcess", # Matric showing inputs and outputs the processes of each plant.
			"fuelPrice", # Prices of fuels (by Time)
			"plantParameters" # Capacities and costs of running plants
		),
		Ident = c(
			NA,
			NA,
			NA,
			NA,
			NA,
			"Op_en4151/fuelPrice",
			NA
		)
	),
	formula = function(...) {
		
		# First, define the optimising function
		optfunction <- function(
			x, # Variables to optimise (activities of power plants).
			amountRequired, # energy demand that must be met by supply (scalar)
			unitprice, # vector of prices per activity unit (by plant).
			unitRequired # vector of required energy produced per activity uniit (by plant).
		) {

			out <- sum(x * unitprice) + (amountRequired - sum(x * unitRequired))^2 * 1E6
			
			return(out)
		}	

		ene <- unkeep(energy, prevresults = TRUE, sources = TRUE) # Modified energy supply and demand
		ene <- oapply(ene, cols = c("Efficiency", "Renovation", "Building"), FUN = sum)
		#requiredEnergy@name <- "requiredEnergy"
		#colnames(requiredEnergy@output)[colnames(requiredEnergy@output) == "energyUseResult"] <- "requiredEnergyResult"
		#colnames(reqen@output)[colnames(reqen@output) == "Heating"] <- "Fuel"
		ene <- oapply(ene * fuelShares, cols = "Heating", FUN = sum)

		if(!is.na(requiredName)) { # Energy balance is only done if requiredName is given

			nondynen <- ene # Non-dynamic energy supply and demand, which is not optimised.
			nondynen@output <- nondynen@output[nondynen@output$Fuel != requiredName , ]

			# Required energy supply and demand, which IS optmised.
			amountRequired <- ene
			if("Fuel" %in% colnames(amountRequired@output)) {
				amountRequired@output <- amountRequired@output[amountRequired@output$Fuel == requiredName , ]
				amountRequired <- unkeep(amountRequired, cols = "Fuel")
			}

			# Operational costs per unit of activity. Fuel price is excluded and calculated elsewhere.
			operationCost <- plantParameters
			operationCost@output <- operationCost@output[operationCost@output$Parameter == "Operation cost" , ]
			operationCost <- unkeep(operationCost, cols = "Parameter")

			# Costs per unit of activity by plant. Fuel price included.
			unitprice <- oapply(energyProcess * fuelPrice * -1, cols = "Fuel", FUN = sum) + operationCost
			
			colnames(unitprice@output)[colnames(unitprice@output) == "plantParametersResult"] <- "unitpriceResult"
			unitprice@output$unitpriceResult <- result(unitprice) # Save result into another column
			
			# Required energy produced per unit activity by plant.
			unitRequired <- energyProcess
			unitRequired@output <- unitRequired@output[unitRequired@output$Fuel == requiredName , ]
			unitRequired <- unkeep(unitRequired, cols = "Fuel")

			# Put all inputs into one ovariable and sort by Plant. Note: important outputs in columns
				# unitpriceResult, requiredEnergyResult, energyProcessResult, and Result.
			temp <- (unitprice * unitRequired * amountRequired * 0 +  
				unkeep(plantParameters, prevresults = TRUE, sources = TRUE)
			)
			indices <- unique(temp@output[temp@marginal & !colnames(temp@output) %in% c("Plant", "Parameter")])

			out <- data.frame()
			
			for(i in 1:nrow(indices)) {
				temp2 <- merge(indices[i , , drop = FALSE], temp@output)
				temp2 <- temp2[order(temp2$Plant) , ]
				temp3 <- optim(
					par			= temp2[temp2$Parameter == "Min" , "Result"], # from plantParameters
					fn 			= optfunction, 
					lower 		= temp2[temp2$Parameter == "Min" , "Result"], # from plantParameters
					upper 		= temp2[temp2$Parameter == "Max" , "Result"], # from plantParameters
					unitprice 	= temp2[temp2$Parameter == "Min" , "unitpriceResult"], # from unitprice
					amountRequired = temp2[temp2$Parameter == "Min" , paste(ene@name, "Result", sep = "")][1], # from amountRequired
					unitRequired = temp2[temp2$Parameter == "Min" , "energyProcessResult"], # from energyProcess
					method = "L-BFGS-B"
				)

				out <- rbind(out, data.frame(
					indices[i , , drop = FALSE],
					Plant = temp2$Plant[temp2$Parameter == "Min"],
					Result = temp3$par
				))
			}
			out <- Ovariable(output = out, marginal = c(rep(TRUE, ncol(out)-1), FALSE))
			out <- out * energyProcess # Optimised fuel usage
#			out <- combine(out, nondynen) # Combine optimised and non-optimised energy demand
		} else {
			out <- ene # If optimisation is not done, take the simple version updated with fuelShares.
		}
		
#		out <- combine(out, nondynsupply) # Combine also non-dynamic energy supply

		return(out)
	}
)

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

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")