mc_simulator.py

import random
import math
import time

NUMBER_OF_SIMULATIONS = 1000

import math
def lognormal_mean(mu, sigma):
    sigma_s = sigma * sigma
    mu_s = mu * mu
    return math.log(mu) - 0.5 * math.log(float(sigma_s)/float(mu_s) + 1.0)
def lognormal_sigma(mu, sigma):
    sigma_s = sigma * sigma
    mu_s = mu * mu
    return math.sqrt(math.log(float(sigma_s)/float(mu_s) + 1.0))
def lognormal_max(mu, sigma):
    mean = lognormal_mean(mu, sigma)
    sigma = lognormal_sigma(mu, sigma)
    return math.e**(mean + 1.96 * sigma)
def lognormal_min(mu, sigma):
    mean = lognormal_mean(mu, sigma)
    sigma = lognormal_sigma(mu, sigma)
    return math.e**(mean - 1.96 * sigma)
parameters_ranges = {
'Cin' : (lognormal_mean(400, 80), lognormal_sigma(400, 80), 'lognormal', lognormal_min(400, 80), lognormal_max(400, 80)),# Values are based upon monitored data
'Co'  : (lognormal_mean(425, 85), lognormal_sigma(425, 85), 'lognormal', lognormal_min(425, 85), lognormal_max(425, 85)),# Values are based upon monitored data
'n'   : (lognormal_mean(7, 4.2), lognormal_sigma(7, 4.2),'lognormal', lognormal_min(7, 4.2), lognormal_max(7, 4.2)),# Values are based upon monitored data
'q'   : (lognormal_mean(2, 1.4), lognormal_sigma(2, 1.4), 'lognormal', lognormal_min(2, 1.4), lognormal_max(2, 1.4)),# Values are based upon monitored data 
'V'   : (lognormal_mean(340, 170), lognormal_sigma(340, 170), 'lognormal', lognormal_min(340, 170), lognormal_max(340, 170)),# Values are based upon monitored data
't'   : (1, 1,  'uniform', 0, 2),
}

parameters_random_values = {}
print parameters_ranges

def conc_equation(Cin, Co, n, q, v , t):
    return ((q/ (n*v)) * (1 - math.e**(-n*t))*1000000) + ((Co - Cin) * math.e**(-n*t)) + Cin

#Where
#Cin = 350-450 ppm
#Co = 350-500 ppm
#n = 2-12 (1/hr)
#q = 0.5-3.5 (m3/hr)
#V  = 100-580 m3
##Cin > Co
for parameter, ranges in parameters_ranges.iteritems():
    if parameter == 'Cin':
        continue
    par1 = ranges[0]
    par2 = ranges[1]
    prob_dist = ranges[2]
    minimum = ranges[3]
    maximum = ranges[4]
    ranges_array = []
    loop_counter = 0
    random.seed(time.time())
    parameters_random_values[parameter] = []
    while (loop_counter < NUMBER_OF_SIMULATIONS):
        if prob_dist == 'uniform':
            random_number = random.uniform(par1, par2)
        if prob_dist == 'normal':
            random_number = random.normalvariate(par1, par2)
        if prob_dist == 'lognormal':
            random_number = random.lognormvariate(par1, par2)
        if minimum > random_number or maximum < random_number:
            continue
        parameters_random_values[parameter].append(random_number)   
        loop_counter = loop_counter + 1

parameter ='Cin'
ranges = parameters_ranges[parameter]
par1 = ranges[0]
par2 = ranges[1]
prob_dist = 'lognormal'
ranges_array = []
loop_counter = 0
minimum = ranges[3]
maximum = ranges[4]

random.seed(time.time())
parameters_random_values[parameter] = []
while (loop_counter < NUMBER_OF_SIMULATIONS):
    if prob_dist == 'uniform':
        upper_limit = parameters_random_values['Co'][loop_counter]
        random_number = random.uniform(par1, par2)
        #Condition on Cin      
        while (random_number >= upper_limit or random_number < 0):
            random_number = random.uniform(par1, par2)
    if prob_dist == 'normal':
        upper_limit = parameters_random_values['Co'][loop_counter]
        random_number = random.normalvariate(par1, par2)
        #Condition on Cin       
        while (random_number >= upper_limit or random_number < 0):
            random_number = random.normalvariate(par1, par2)
    if prob_dist == 'lognormal':
        upper_limit = parameters_random_values['Co'][loop_counter]
        random_number = random.lognormvariate(par1, par2)
        #Condition on Cin     
        while (random_number >= upper_limit or random_number < 0):
            random_number = random.lognormvariate(par1, par2)
            
    if minimum > random_number or maximum < random_number:
        continue
    parameters_random_values[parameter].append(random_number)   
    loop_counter = loop_counter + 1

loop_counter = 0
output_file = open('C:/Python27/conc_monte_carlo_simulation.csv', 'w')
output_file.write("Cin, Co, n, q, V, T, SIMULATION RESULT(C)\n")

while (loop_counter < NUMBER_OF_SIMULATIONS):  
    Cin = parameters_random_values['Cin'][loop_counter]
    Co  = parameters_random_values['Co'][loop_counter]
    n   = parameters_random_values['n'][loop_counter]
    q   = parameters_random_values['q'][loop_counter]
    V   = parameters_random_values['V'][loop_counter]
    t   = parameters_random_values['t'][loop_counter]
    output_file.write(",".join([str(Cin), str(Co), str(n), str(q), str(V) , str(t), str(conc_equation(Cin, Co, n, q, V , t))]) + "\n")
    loop_counter = loop_counter + 1

output_file.close()

##In this part of the programme, the output excel file is taken as input and one more excel file with cout (i.e. output concentration) values and frequency is generated.

result_file = open('C:/Python27/conc_monte_carlo_simulation.csv')

results = result_file.read()
result_list = results.split('\n')

del(result_list[0])
del(result_list[NUMBER_OF_SIMULATIONS])

result_list_list = []
for result in result_list:
    result_list_list.append(result.split(','))
    cout = result_list_list[-1][-1]
    result_list_list[-1].append(math.floor(float(cout)))

bucketed_results = {}
for result in result_list_list:
    cout = result[-1]
    if cout not in bucketed_results:
        bucketed_results[cout]=0
    bucketed_results[cout] = bucketed_results[cout] + 1

output_histogram = open('C:/Python27/output_histogram.csv', 'w')

for cout_value, frequency in bucketed_results.iteritems():
    output_histogram.write(str(cout_value) + ',' + str(frequency) + '\n')
   
output_histogram.close()