# -*- coding: utf-8 -*-
"""
@author: Vaclav Miklas

Nazev: TECHNICKO-EKONOMICKA ANALYZA
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import bisect
import pandas as pd

# VYSLEDKY Z MODELU PRO VZOROVOU BPS O EL. VYKONU CCA 400 kWe:
mFug = 3041.54 # ... hmotnostni tok [kg/h]
susFug = 0.03 # ... susina ve fugatu [–]
# Pri 20 % vyuziti tepla na vytapeni fermentoru produkce "odsolene" vody priblizne
# odpovida 50 m3 – investicni naklady na MSF o teto kapacite (vyrobce PBS)
# dle dizertacni prace Dr. Hajka jsou 45000 €, coz aktualne odpovida 1,15 mil. Kc.
# Tato castka byla s ohledem na "spinavou" aplikaci navysena o 20 %.
C = [1725000, 0, 0] # ... investicni naklady [Kc]
A = [199.93, 174.79, 149.73] # ... celkova teplosmenna plocha MSF [m2]
E = [24.55, 21.69, 18.91] # ... spotreba elektricke energie [kWe]
mConc = [940.96, 1204.41, 1467.35] # ... hmotnostni tok koncentratu [kg/h]

# DALSI PARAMETRY:
r = 0.10 # [–] ... pro stredne rizikovy projekt je zvolena nominalni diskontni 
        #         sazba 10 %
m = 0.025 # [–] ... ocekavana mira inflace 2,5 %
#vykup = 3.55 # [Kc/kWh] ... vykupni cena za elektrinu (odpovida napr. BPS do
             #              550 kWe uvedene do provozu ve 2013)
rk = ((1+r)/(1+m)-1) # [–] realna (korigovana) diskontni sazba
hod = 8200 # [h] ... rocni fond provozni doby

#==============================================================================
# DOPOCET INVESTICNICH NAKLADU:
# Dle Product & Process Design Principles je vztah pro vetsinu klasickych odparek
# ve tvaru C = k*A^0,53. MSF zde uvedena neni, avsak obecne pro zarizeni,
# u nichz inv. naklady zavisi na teplosmenne plose, lezi typicky exponent v
# rozmezi 0,4 a 0,6. Jako referencni je vzata cena viz vyse:

q = C[0]/A[0]**0.53
for i in (1,2):
    C[i] = q*A[i]**0.53
    
#==============================================================================
# VYPOCET NAVRATNOSTI Z DISKONTOVANEHO KUMULATIVNIHO CASHFLOW:

def payback(doprava, vykup, C, A, E, mConc):
    # Zavedeni promennych:
    uspora = [0]*len(C) # uspora na doprave
    ztrat_el = [0]*len(C) # ztratova elektrina, ktera by byla k prodeji
    CF = [0]*len(C) # cashflow
    DCF = np.zeros((3,100)) # diskontovane cashflow
    DCCF = np.zeros((3,100)) # diskontovane kumulativni cashflow
    PB = [0]*3 # navratnost investice
    # Pro kazdy ze tri pripadu:
    for i in range(len(C)):
        # Se stanovi uspora, ztratova elektrina a cashflow:
        uspora[i] = doprava*(mFug-mConc[i])*hod/1000
        ztrat_el[i] = E[i]*vykup*hod
        CF[i] = uspora[i]-ztrat_el[i]
        # Predpoklada se, ze tyto hodnoty se mezirocne nemeni.
        for j in range (1,51,1):
            # Pro kazdy rok (nejvyse 100) je napocitane diskontovane CF a
            # diskontovane kumulativni CF.
            DCF[i][j] = CF[i]/(1+rk)**j
            DCCF[i][j] = DCCF[i][j-1]+DCF[i][j]
            # Jakmile DCCF je vetsi nez investice...
            if DCCF[i][j]/C[i] > 1:
                # ...tak je zaznemana doba navratnosti a prerusen cyklus:
                PB[i] = j-(DCCF[i][j]-C[i])/(DCCF[i][j]-DCCF[i][j-1])
                break
                # Protoze DCCF neni spojita velicina, pro casti roku je
                # pouzita linearni aproximace.

    return PB
        
def payback_iter(doprava):
    # Zavedeni promennych:
    uspora = [0]*len(C) # uspora na doprave
    ztrat_el = [0]*len(C) # ztratova elektrina, ktera by byla k prodeji
    CF = [0]*len(C) # cashflow
    DCF = np.zeros((3,100)) # diskontovane cashflow
    DCCF = np.zeros((3,100)) # diskontovane kumulativni cashflow
    PB = [0]*3 # navratnost investice
    # Pro dany pripad se stanovi uspora, ztratova elektrina a cashflow:
    uspora[i] = doprava*(mFug-mConc[i])*hod/1000
    ztrat_el[i] = E[i]*vykup*hod
    CF[i] = uspora[i]-ztrat_el[i]
    # Predpoklada se, ze tyto hodnoty se mezirocne nemeni.
    for j in range (1,51,1):
        # Pro kazdy rok (nejvyse 100) je napocitane diskontovane CF a
        # diskontovane kumulativni CF.
        DCF[i][j] = CF[i]/(1+rk)**j
        DCCF[i][j] = DCCF[i][j-1]+DCF[i][j]
        # Jakmile DCCF je vetsi nez investice...
        if DCCF[i][j]/C[i] > 1:
            # ...tak je zaznemana doba navratnosti a prerusen cyklus:
            PB[i] = j-(DCCF[i][j]-C[i])/(DCCF[i][j]-DCCF[i][j-1])
            break
            # Protoze DCCF neni spojita velicina, pro casti roku je
            # pouzita linearni aproximace.
    return (PB[i]-poz_navr)

writer = pd.ExcelWriter('Navratnosti.xlsx', engine='xlsxwriter')
df = {}
for poz_navr in (2,5,10):
    cena = [0]*3
    k = 0
    for i in range(3):
        pomocna = []
        for vykup in (0.9,3.55,4.12):
            k += 1
            if k==1:
                low = 23
            elif k==2:
                low = 55
            elif k==3:
                low = 62
            pomocna.append(bisect(payback_iter,low,500))
        k = 0
        cena[i]=pomocna
    df[poz_navr] = pd.DataFrame(data=cena, index = (80,70,60), columns=(0.9,3.55,4.12))
    df[poz_navr].to_excel(writer, sheet_name='Navratnost '+str(poz_navr))
writer.save()
    

#==============================================================================
# VYKRESLENI GRAFU:    
pom = []
navratnosti = [0]*len(C)

k = 0

for vykup in (0.9,3.55,4.12):
    a = True
    x_vec = []
    k += 1
    if k==1:
        low = 23
        high = 91
    elif k==2:
        low = 55
        high = 121
    elif k==3:
        low = 62
        high = 121
    for i in range(len(C)):
        pom = []
        for doprava in range (low,high,1):    
            pom.append(payback(doprava,vykup,C,A,E,mConc)[i])
            if a == True:
                x_vec.append(doprava)
            navratnosti[i] = pom
        a = False
            
    plt.figure ()
    plt.rc('font', family='serif')
    
    plt.plot(x_vec, navratnosti[0], 'k-', label=u'80 % vyrobeného tepla k dispozici')
    plt.plot(x_vec, navratnosti[1], 'k-.', label=u'70 % vyrobeného tepla k dispozici')
    plt.plot(x_vec, navratnosti[2], 'k--', label=u'60 % vyrobeného tepla k dispozici')
   
    plt.title(u'Výkupní cena elektřiny '+str(vykup).replace('.',',')+' Kč/kWh')
    plt.xlabel(u'Cena za dopravu [Kč/t]')
    plt.ylabel(u'Návratnost [rok]')
    plt.legend()
    # plt.show()
    plt.savefig(str(vykup).replace('.','_')+'.svg', dpi=1200, transparent=True)
    