Моделирование Монте-Карло с треугольным и нормальным распределением плотности вероятности с использованием NumPy - PullRequest
0 голосов
/ 10 сентября 2018

Я пытаюсь выполнить моделирование по методу Монте-Карло, чтобы вычислить неопределенность в затратах на электроэнергию для системы с тепловым насосом.У меня есть несколько входных параметров (COP, затраты на электроэнергию), которые имеют треугольное распределение вероятностей.Общие затраты на электроэнергию состоят из суммы рассчитанных затрат трех подкомпонентов (тепловых насосов и насосов) и имеют (приблизительно) нормальное распределение вероятностей.

Мне было интересно, правильно ли я выполняю симуляцию MC,Поскольку мне приходится зацикливать эту симуляцию MC на более чем 70 различных системах тепловых насосов, мне также интересно, есть ли более быстрый метод.

Поскольку я абсолютный приверженец кодирования, прошу прощения за мой грязный код.

Я благодарен за любую помощь!

Мой код:

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)

1 Ответ

0 голосов
/ 10 сентября 2018

Отлично, у вас есть прототип для вашего кода!

Некоторые впечатления о структуре кода и удобочитаемости:

  • самое быстрое улучшение - это разделение функций и части сценариев, это позволяет разбить ваш алгоритм на простые тестируемые блоки и сохранить контроль вычислений в одном месте с последующим построением графика

  • некоторые повторяющиеся коды могут переходить к собственным функциям

  • это окупается, чтобы придерживаться более широко принятого соглашения об именах (PEP8), таким образом, люди меньше удивляются стилю и могут уделять больше внимания содержанию кода. В частности, вы обычно называете функции строчными буквами подчеркивания def do_somtheing():, а UPPERCASE зарезервировано для констант.

Необходимо иметь возможность запускать свой код для проверки ускорений, см. Комментарии выше о том, чего не хватает.

Обновление более полного кода, о котором идет речь, некоторые дополнительные комментарии:

  • пусть функции являются функциями, не передавайте аргументы функции как глобальные
  • рассмотрите возможность отделения структуры вашей модели (ваших уравнений) от параметров (фиксированных входных данных и начальных значений), в более долгосрочной перспективе она окупится
  • избегайте «магических чисел», помещайте их в константы (например, 866) или передавайте явно как аргументы.

Вот рефракция для рассмотрения:

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)  

Следующим неотложным намерением является написание HeatPump класса, но я предлагаю придерживаться до тех пор, пока вы можете - это обычно заставляет нас думать о классе состояние и методы лучше.

Мне также кажется, что значения параметров для параметров могут быть не независимыми, в некоторых экспериментах вы можете извлечь их из совместного распределения.

...