Monte Carlo Simulation with triangular and normal probability density distibution using numpy

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)


triangular()

sim

minval = np.min(Eout[np.nonzero(sim)])

1 Answer
1

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 fucntions

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_somtheing(): and UPPERCASE is reserved for constants.

def do_somtheing():

UPPERCASE

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

Here is a refctoring 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 directrly, this is a very readable way
# make a vector or random arguements
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 fucntions as long as you can - that usually makes our thinking about a class
state and methods better.

HeatPump

I also appears the seeded values for parameters may not be independent, in some experiment you can draw them from a joint distribution.

import *

Thanks for contributing an answer to Stack Overflow!

But avoid …

To learn more, see our tips on writing great answers.

Required, but never shown

Required, but never shown

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.