I am trying to perform a Monte Carlo Simulation in order to calculate the uncertainty in the electricity costs of a heat pump system. I have several input parameters (COPs, electricity costs), which are of a triangular probability distribution. The total electricity costs are composed of the sum of the calculated costs of the three subcomponents (heatpump and pumps) and are of an (approximately) normal probability distribution.
I was wondering if I am performing the MC simulation correctly. Since I have to loop this MC simulation over 70 different heat pump systems, I am also wondering if there is a faster method.
As I am an absolute greenhorn in coding, please apologize for my messy code.
I am thankful for any help!
My code:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import triangular
N = 1_000_000
def energy_output(coef_performance, energy_input):
return energy_input * coef_performance / (coef_performance - 1)
COP_DISTRIBUTION_PARAM = dict(left=4, mode=4.5, right=5)
def seed_cop():
return triangular(**COP_DISTRIBUTION_PARAM )
INPUT_ENERGY_HEATING = 866
INPUT_ENERGY_COOLING = 912
def random_energy_output():
return energy_output(seed_cop(), energy_input=INPUT_ENERGY_HEATING)
energy_outputs = [random_energy_output() for _ in range(N)]
a = min(energy_outputs)
b = max(energy_outputs)
med = np.median(energy_outputs)
############################
def elec_costs_heatpump(elec_costs, coef_performance,energy_output):
return energy_output * 1000 / coef_performance * elec_costs
ELEC_DISTRIBUTION_PARAM = dict(left=0.14, mode=0.15, right=0.16)
def seed_elec():
return triangular(**ELEC_DISTRIBUTION_PARAM )
HP_OUTPUT_DISTRIBUTION_PARAM = dict(left=a, mode=med, right=b)
def seed_output():
return triangular(**HP_OUTPUT_DISTRIBUTION_PARAM )
def random_elec_costs_heatpump():
return elec_costs_heatpump(seed_elec(),seed_cop(),seed_output() )
elec_costs_heatpump = [random_elec_costs_heatpump() for _ in range(N)]
mean_hp = np.mean(elec_costs_heatpump)
std_hp = np.std(elec_costs_heatpump)
############################
def elec_costs_coldpump(elec_costs, coef_performance_pump,energy_input):
return energy_input * 1000 / coef_performance_pump * elec_costs
COP_PUMP_DISTRIBUTION_PARAM = dict(left=35, mode=40, right=45)
def seed_cop_pump():
return triangular(**COP_PUMP_DISTRIBUTION_PARAM )
def random_elec_costs_coldpump():
return elec_costs_coldpump(seed_elec(),seed_cop_pump(), energy_input=INPUT_ENERGY_COOLING)
elec_costs_coldpump = [random_elec_costs_coldpump() for _ in range(N)]
mean_cp = np.mean(elec_costs_coldpump)
sdt_cp = np.std(elec_costs_coldpump)
#########################
def elec_costs_warmpump(elec_costs, coef_performance_pump,energy_input):
return energy_input * 1000 / coef_performance_pump * elec_costs
def random_elec_costs_warmpump():
return elec_costs_warmpump(seed_elec(),seed_cop_pump(), energy_input=INPUT_ENERGY_HEATING)
elec_costs_warmpump = [random_elec_costs_warmpump() for _ in range(N)]
mean_wp = np.mean(elec_costs_warmpump)
sdt_wp = np.std(elec_costs_warmpump)
#########################
def total_costs(costs_heatpump, costs_coldpump, costs_warmpump):
return costs_heatpump + costs_coldpump + costs_warmpump
ELEC_COSTS_HEATPUMP_PARAM = dict(loc=mean_hp, scale=sdt_hp)
def seed_costs_hp():
return np.random.normal(**ELEC_COSTS_HEATPUMP_PARAM )
ELEC_COSTS_COLDPUMP_PARAM = dict(loc=mean_cp, scale=sdt_cp)
def seed_costs_cp():
return np.random.normal(**ELEC_COSTS_COLDPUMP_PARAM )
ELEC_COSTS_WARMPUMP_PARAM = dict(loc=mean_wp,scale=sdt_wp)
def seed_cost_wp():
return np.random.normal(**ELEC_COSTS_WARMPUMP_PARAM )
def random_total_costs():
return seed_costs_hp(), seed_costs_cp(), seed_cost_wp()
total_costs = [random_total_costs() for _ in range(N)]
print(total_costs)
#Plot = plt.hist(total_costs, bins=75, density=True)
Great you have a prototype for your code!
Some impressions on code structure and readability:
the quickest improvement is separating functions and scripting part, this allows to split your algorithm into simple testable blocks and keep control of calculations in one place, followed by plotting
some repeated code can go to own functions
it pays back to stick to more widely accepted naming convention (PEP8), this way people are less astonished less with style and can devote more attention to code substance. Specifically, you normally name functions lowercase underscore def do_something(): and UPPERCASE is reserved for constants.
Need to be able to run your code to check for speedups, see comments above about what is lacking.
Update on more complete code in question, some additional comments:
866) or pass explicitly as args.Here is a refactoring to consider:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import triangular
N = 1_000_000
#def EnergyHP():
# COP_Hp = triangular(4.0, 4.5, 5) # COP of heatpump
# Ein = 866 # Energy input heatpump
# return Ein * COP_Hp /(COP_Hp-1) # Calculation of Energy output of heatpump
#
#Eout = np.zeros(N)
#for i in range(N):
# Eout[i] = EnergyHP()
#minval = np.min(Eout[np.nonzero(Eout)])
#maxval = np.max(Eout[np.nonzero(Eout)])
#Medi= np.median(Eout, axis=0)
def energy_output(coef_performance, energy_input):
"""Return energy output of heatpump given efficiency and input.
Args:
coef_performance - <description, units of measurement>
energy_input - <description, units of measurement>
"""
return energy_input * coef_performance / (coef_performance - 1)
# you will use this variable again, so we put it into a function to recylce
COP_DISTRIBUTION_PARAM = dict(left=4, mode=4.5, right=5)
def seed_cop():
"""Get random value for coef of performance."""
return triangular(**COP_DISTRIBUTION_PARAM )
# rename it if it is a constant, pass as an argument is it is not
SOME_MEANINGFUL_CONSTANT = 866
def random_energy_output():
return energy_output(seed_cop(), energy_input=SOME_MEANINGFUL_CONSTANT)
# pure python list and metrics
energy_outputs = [random_energy_output() for _ in range(N)]
a = min(energy_outputs)
b = max(energy_outputs)
med = np.median(energy_outputs)
# above does this does not use numpy, but you can convert it to vector
energy_outputs_np = np.array(energy_outputs)
# or you can construct np array directly, this is a very readable way
# make a vector or random arguments
cop_seeded = triangular(**COP_DISTRIBUTION_PARAM, size=N)
# apply function
energy_outputs_np_2 = energy_output(cop_seeded, SOME_MEANINGFUL_CONSTANT)
The next pressing intent is writing a HeatPump class, but I suggest sticking
to functions as long as you can - that usually makes our thinking about a class
state and methods better.
I also appears the seeded values for parameters may not be independent, in some experiment you can draw them from a joint distribution.