# -*- coding: utf-8 -*-
"""
Created on Mon Sep 28 2020

"""

import openpyxl
import gurobipy as gp
from gurobipy import GRB

### define which optimizations are to be done
co2       = True #set True for co2 optimization
econ_lin  = True #set True for economic optimization with linear cost structure
econ_nlin = True #set True for economic optimization with non-linear cost structure

#---------------------------------------------------------------------------

### read excel workbook
par_wb    = openpyxl.load_workbook("input parameters.xlsx", data_only=True)
par_sheet = par_wb["parameters"]


### create parameters from excel sheet
IR        = float(par_sheet["F3"].value)/100     #) Discount rate
P_PU      = float(par_sheet["F4"].value)         #) Pulp - Annual production 
P_EL      = float(par_sheet["F5"].value)         #) Electric power – Annual excess production

D_O       = float(par_sheet["F6"].value)         #) Oxygen – Annual demand
PP_O      = float(par_sheet["F7"].value)         #) Oxygen – Purchasing price

D_T       = float(par_sheet["F8"].value)         #) Trucks for wood supply – Needed number of trucks
D_TKM     = float(par_sheet["F9"].value)         #) Trucks for wood supply – Annual driven distance
FC_T      = float(par_sheet["F10"].value)        #) Trucks for wood supply – Fuel consumption per truck
PP_FT     = float(par_sheet["F11"].value)        #) Trucks for wood supply – Fuel price
EM_T      = float(par_sheet["F12"].value)        #) Trucks for wood supply – CO2 emissions per truck

D_F       = float(par_sheet["F13"].value)        #) Forklifts – Needed number of forklifts
D_FKM     = float(par_sheet["F14"].value)        #) Forklifts – Annual driven distance
FC_F      = float(par_sheet["F15"].value)        #) Forklifts – Fuel consumption per forklift
PP_FF     = float(par_sheet["F16"].value)        #) Forklifts – Fuel price
EM_F      = float(par_sheet["F17"].value)        #) Forklifts – CO2 emissions per forklift

D_NGLK    = float(par_sheet["F18"].value)        #) Lime kiln – Annual natural gas consumption
PP_NG     = float(par_sheet["F19"].value)        #) Lime kiln – Price of natural gas
EM_LK     = float(par_sheet["F20"].value)        #) Lime kiln – Annual CO2 emissions
R_HLK     = float(par_sheet["F21"].value)/100    #) Lime kiln – Max. hydrogen injection rate

SP_O      = float(par_sheet["F22"].value)        #) Oxygen – Sale potential
PS_O      = float(par_sheet["F23"].value)        #) Oxygen – Sale price

SP_HC     = float(par_sheet["F24"].value)        #) Hydrogen – Cylinder sale potential
PS_HC     = float(par_sheet["F25"].value)        #) Hydrogen – Cylinder sale price
SP_HI     = float(par_sheet["F26"].value)        #) Hydrogen – Gas grid injection sale potential
PS_HI     = float(par_sheet["F27"].value)        #) Hydrogen – Gas grid injection sale price
CC_HT     = float(par_sheet["F28"].value)        #) Hydrogen – Replacement truck price
LT_HT     = float(par_sheet["F29"].value)        #) Hydrogen – Replacement truck lifetime
FC_HT     = float(par_sheet["F30"].value)        #) Hydrogen – Replacement truck hydrogen consumption
CC_HF     = float(par_sheet["F31"].value)        #) Hydrogen – Replacement forklift price
LT_HF     = float(par_sheet["F32"].value)        #) Hydrogen – Replacement forklift lifetime
FC_HF     = float(par_sheet["F33"].value)        #) Hydrogen – Replacement forklift hydrogen consumption
CC_HFI    = float(par_sheet["F34"].value)        #) Hydrogen – Vehicle refilling station capital cost
CO_HFI    = float(par_sheet["F35"].value)/100    #) Hydrogen – Vehicle refilling station annual operating cost
LT_HFI    = float(par_sheet["F36"].value)        #) Hydrogen – Vehicle refilling station lifetime
CC_KHFI   = float(par_sheet["F37"].value)        #) Hydrogen – Known capital cost of known filling station with capacity CA_KHFI
CA_KHFI   = float(par_sheet["F38"].value)        #) Hydrogen – Capacity of known filling station
N_KHFI    = float(par_sheet["F39"].value)        #) Hydrogen – Filling station’s exponent for six-tenth rule

SP_CO2    = float(par_sheet["F40"].value)        #) CO2 – Sale potential
PS_CO2    = float(par_sheet["F41"].value)        #) CO2 – Sale price

SP_LI     = float(par_sheet["F42"].value)        #) Lignin – Sale potential
PS_LI     = float(par_sheet["F43"].value)        #) Lignin – Sale price

SP_MEC    = float(par_sheet["F44"].value)        #) Methane – Sale potential
PS_MEC    = float(par_sheet["F45"].value)        #) Methane – Sale price
SP_MEI    = float(par_sheet["F46"].value)        #) Methane – Grid injection potential
PS_MEI    = float(par_sheet["F47"].value)        #) Methane – Grid injection sale price

SP_ML     = float(par_sheet["F48"].value)        #) Methanol – Sale potential
PS_ML     = float(par_sheet["F49"].value)        #) Methanol – Sale price

PR_OPE    = float(par_sheet["F50"].value)        #) Electrolysis – Oxygen produced per kWh electric power
PR_HPE    = float(par_sheet["F51"].value)        #) Electrolysis – Hydrogen produced per kWh electric power
CC_EL     = float(par_sheet["F52"].value)        #) Electrolysis – Capital cost
CO_EL     = float(par_sheet["F53"].value)/100    #) Electrolysis – Annual operating cost
LT_EL     = float(par_sheet["F54"].value)        #) Electrolysis – Lifetime
OH_EL     = float(par_sheet["F55"].value)        #) Electrolysis – Annual operating hours
CC_SEL    = float(par_sheet["F56"].value)        #) Electrolysis – Stack replacement cost
LT_SEL    = float(par_sheet["F57"].value)        #) Electrolysis – Stack lifetime
CC_KEL    = float(par_sheet["F58"].value)/100    #) Electrolysis – Known capital cost of known plant with capacity CA_KEL
CA_KEL    = float(par_sheet["F59"].value)        #) Electrolysis – Capacity of known plant
N_KEL     = float(par_sheet["F60"].value)        #) Electrolysis – Exponent for six-tenth rule

PR_MEPCO2 = float(par_sheet["F61"].value)        #) Methanation – Methane produced per t CO2 
PR_MEPH   = float(par_sheet["F62"].value)        #) Methanation – Methane produced per kg hydrogen
CC_ME     = float(par_sheet["F63"].value)        #) Methanation – Capital cost
CO_ME     = float(par_sheet["F64"].value)/100    #) Methanation – Annual operating cost
LT_ME     = float(par_sheet["F65"].value)        #) Methanation – Lifetime
OH_ME     = float(par_sheet["F66"].value)        #) Methanation – Annual operating hours
CC_KME    = float(par_sheet["F67"].value)        #) Methanation – Known capital cost of known plant with capacity CA_KME
CA_KME    = float(par_sheet["F68"].value)        #) Methanation – Capacity of known plant
N_KME     = float(par_sheet["F69"].value)        #) Methanation – Exponent for six-tenths rule

PR_MLPCO2 = float(par_sheet["F70"].value)        #) Methanol synthesis – Methanol produced per t CO2 
PR_MLPH   = float(par_sheet["F71"].value)        #) Methanol synthesis – Methanol produced per kg hydrogen
CC_ML     = float(par_sheet["F72"].value)        #) Methanol synthesis – Capital cost
CO_ML     = float(par_sheet["F73"].value)/100    #) Methanol synthesis – Annual operating cost
LT_ML     = float(par_sheet["F74"].value)        #) Methanol synthesis – Lifetime
OH_ML     = float(par_sheet["F75"].value)        #) Methanol synthesis – Annual operating hours
CC_KML    = float(par_sheet["F76"].value)        #) Methanol synthesis – Known capital cost of known plant with capacity CA_KML
CA_KML    = float(par_sheet["F77"].value)        #) Methanol synthesis – Capacity of known plant
N_KML     = float(par_sheet["F78"].value)        #) Methanol synthesis – Exponent for six-tenths rule

CC_CC     = float(par_sheet["F79"].value)        #) CCS – Capital cost
CO_CC     = float(par_sheet["F80"].value)/100    #) CCS – Annual operating cost
LT_CC     = float(par_sheet["F81"].value)        #) CCS – Lifetime
CC_KCC    = float(par_sheet["F82"].value)        #) CCS – Known capital cost of known plant with capacity CAKCC
CA_KCC    = float(par_sheet["F83"].value)        #) CCS – Capacity of known plant
N_KCC     = float(par_sheet["F84"].value)        #) CCS – Exponent for six-tenth rule

PR_LIPCO2 = float(par_sheet["F85"].value)        #) Lignoboost – Lignin produced per t CO2
CC_LI     = float(par_sheet["F86"].value)        #) Lignoboost – Capital cost
CO_LI     = float(par_sheet["F87"].value)/100    #) Lignoboost – Annual operating cost
LT_LI     = float(par_sheet["F88"].value)        #) Lignoboost – Lifetime
CC_KLI    = float(par_sheet["F89"].value)        #) Lignoboost – Known capital cost of known plant with capacity CA_KLI
CA_KLI    = float(par_sheet["F90"].value)        #) Lignoboost – Capacity of known plant
N_KLI     = float(par_sheet["F91"].value)        #) Lignoboost – Exponent for six-tenth rule




### calculate annuities
A_HT  = (1 - 1 /(1 + IR)**LT_HT) /IR     #hydrogen truck
A_HF  = (1 - 1 /(1 + IR)**LT_HF) /IR     #hydrogen forklift
A_HFI = (1 - 1 /(1 + IR)**LT_HFI) /IR    #hydrogen filling station
A_EL  = (1 - 1 /(1 + IR)**LT_EL) /IR     #electrolysis system
A_SEL = (1 - 1 /(1 + IR)**LT_SEL) /IR    #electrolysis stack
A_CC  = (1 - 1 /(1 + IR)**LT_CC) /IR     #CCS system
A_LI  = (1 - 1 /(1 + IR)**LT_LI) /IR     #lignoboost system
A_ME  = (1 - 1 /(1 + IR)**LT_ME) /IR     #methanation system
A_ML  = (1 - 1 /(1 + IR)**LT_ML) /IR     #methanol synthesis system

#---------------------------------------------------------------------------

### optimize function, optimizes the optimizations that are set as True
def opt(co2, econ_lin, econ_nlin):
    
    #create solution workbook and copy input parameter sheet to it
    sol_wb = par_wb
    
    #execute co2 optimization if wanted and write result to solution workbook
    if co2 == True:
        co2_sol   = co2_opt()
        co2_sheet = "CO2 emission optimization"
        write(sol_wb=sol_wb, sheetname=co2_sheet, opt_sol=co2_sol)
        sol_wb[co2_sheet]["C3"] = co2_sol.objVal
        sol_wb[co2_sheet]["C4"] = econ_lin_fct(co2_sol)
        sol_wb[co2_sheet]["C5"] = econ_nlin_fct(co2_sol)
        
    #execute linear economic optimization if wanted and write result to solution workbook
    if econ_lin == True:
        econ_lin_sol   = econ_lin_opt()
        econ_lin_sheet = "Linear economic optimization"
        write(sol_wb=sol_wb, sheetname=econ_lin_sheet, opt_sol=econ_lin_sol)
        sol_wb[econ_lin_sheet]["C3"] = co2_fct(econ_lin_sol)
        sol_wb[econ_lin_sheet]["C4"] = econ_lin_sol.objVal
        sol_wb[econ_lin_sheet]["C5"] = econ_nlin_fct(econ_lin_sol)
    
    #execute non-linear economic optimization if wanted and write result to solution workbook
    if econ_nlin == True:
        econ_nlin_sol   = econ_nlin_opt()
        econ_nlin_sheet = "Non-linear economic optimization"
        write(sol_wb=sol_wb, sheetname="Non-linear economic optimization", opt_sol=econ_nlin_sol)
        sol_wb[econ_nlin_sheet]["C3"] = co2_fct(econ_nlin_sol)
        sol_wb[econ_nlin_sheet]["C4"] = econ_lin_fct(econ_nlin_sol)
        sol_wb[econ_nlin_sheet]["C5"] = econ_nlin_sol.objVal
    
    sol_wb.save("solution.xlsx")
    
    return True

#writer function, writes results of optimization to workbook
def write(sol_wb, sheetname, opt_sol):
    opt_sol.getVars()
    write_sheet       = sol_wb.create_sheet(title=sheetname)
    
    # Set column width
    write_sheet.column_dimensions['A'].width = 60.0
    write_sheet.column_dimensions['B'].width = 15.0
    write_sheet.column_dimensions['C'].width = 15.0
    
    write_sheet["A1"] = sheetname
    write_sheet["A2"] = "Description"
    write_sheet["B2"] = "Unit"
    write_sheet["C2"] = "Value"
    write_sheet["A3"] = "Annual CO2 emissions saved"
    write_sheet["B3"] = "t/y"
    write_sheet["A4"] = "Annual economic result with linear cost structure"
    write_sheet["B4"] = "€/y"
    write_sheet["A5"] = "Annual economic result with non-linear cost structure"
    write_sheet["B5"] = "€/y"    
    
    units = ["t/y", "t/y", "kg/y", "kg/y", "No. of trucks", "No. of forklifts", "kWh/y", "kg/y", "t/y", "t/y", "t/y", "kWh/y","kWh/y","kWh/y", "kWh/y"]
    i=6
    j=0
    for v in opt_sol.getVars():
        write_sheet["A"+str(i)] = v.varName
        write_sheet["B"+str(i)] = units[j]
        write_sheet["C"+str(i)] = v.x
        i=i+1
        j=j+1
        
        #exit loop before helper variables from economic optimization are included
        if j >= 15:
            return True    
    return True
    
#co2 optimization
def co2_opt():
    print("-----CO2 emission optimization-----")
    
    #create a new model
    m = gp.Model("co2")
    
    #create decision variables
    o_s   = m.addVar(name = "Oxygen - Annual sold amount")        # Oxygen - Annual sold amount
    o_i   = m.addVar(name = "Oxygen - Annual internally used amount")        # Oxygen - Annual internally used amount
    h_s   = m.addVar(name = "Hydrogen - Annual sold amount")        # Hydrogen - Annual sold amount
    h_ig  = m.addVar(name = "Hydrogen - Annual amount injected in gas grid")        # Hydrogen - Annual amount injected in gas grid
    h_t   = m.addVar(name = "Hydrogen - Number of trucks replaced with fuel cell vehicles", vtype = GRB.INTEGER)        # Hydrogen - Number of trucks replaced with fuel cell vehicles
    h_f   = m.addVar(name = "Hydrogen - Number of forklifts replaced with fuel cell vehicles", vtype = GRB.INTEGER)        # Hydrogen - Number of forklifts replaced with fuel cell vehicles
    h_il  = m.addVar(name = "Hydrogen - Annual amount injected in lime kiln")        # Hydrogen - Annual amount injected in lime kiln
    h_to  = m.addVar(name = "Hydrogen - Total amount needed")        # Hydrogen - Total amount needed
    c_s   = m.addVar(name = "CO2 - Annual sold amount")        # CO2 - Annual sold amount
    c_to  = m.addVar(name = "CO2 - Total amount needed")        # CO2 - Total amount needed
    l_s   = m.addVar(name = "Lignin - Annual sold amount")        # Lignin - Annual sold amount
    m_es  = m.addVar(name = "Methane - Annual sold amount")        # Methane - Annual sold amount
    m_eig = m.addVar(name = "Methane - Annual amount injected in gas grid")        # Methane - Annual amount injected in gas grid
    m_eil = m.addVar(name = "Methane - Annual amount injected in lime kiln")        # Methane - Annual amount injected in lime kiln
    m_ls  = m.addVar(name = "Methanol - Annual sold amount")        # Methanol - Annual sold amount
    
    #set objective function
    #obj = EM_T * h_t /1000000
    #obj += EM_F * h_f /1000000
    obj =  h_t * EM_T * D_TKM / 1000000 / D_T
    obj += h_f * EM_F * D_FKM / 1000000 / D_F
    obj += EM_LK / D_NGLK * h_il
    obj += c_s
    m.setObjective(obj, GRB.MAXIMIZE)
    
    #add constraints
    m.addConstr(h_t *FC_HT *D_TKM /100 /D_T + h_f * FC_HF *D_FKM /100 /D_F + h_il /39.45 <= P_EL * 1000 * PR_HPE * 0.08)
    m.addConstr(h_t  <= D_T)
    m.addConstr(h_f  <= D_F)
    m.addConstr(h_il <= D_NGLK * R_HLK)
    m.addConstr(c_s  <= EM_LK)
    m.addConstr(c_s  <= SP_CO2)
    
    #optimize model
    m.optimize()
    
    for v in m.getVars():
        print('%s %g' % (v.varName, v.x))

    print('Obj: %g' % m.objVal)
        
    return m

#economic optimization with linear cost structure
def econ_lin_opt():
    print("-----Economic optimization with linear cost structure-----")
    
    #create a new model
    m = gp.Model("co2")
    
    #create decision variables
    o_s   = m.addVar(name = "Oxygen - Annual sold amount")        # Oxygen - Annual sold amount
    o_i   = m.addVar(name = "Oxygen - Annual internally used amount")        # Oxygen - Annual internally used amount
    h_s   = m.addVar(name = "Hydrogen - Annual sold amount")        # Hydrogen - Annual sold amount
    h_ig  = m.addVar(name = "Hydrogen - Annual amount injected in gas grid")        # Hydrogen - Annual amount injected in gas grid
    h_t   = m.addVar(name = "Hydrogen - Number of trucks replaced with fuel cell vehicles", vtype = GRB.INTEGER)        # Hydrogen - Number of trucks replaced with fuel cell vehicles
    h_f   = m.addVar(name = "Hydrogen - Number of forklifts replaced with fuel cell vehicles", vtype = GRB.INTEGER)        # Hydrogen - Number of forklifts replaced with fuel cell vehicles
    h_il  = m.addVar(name = "Hydrogen - Annual amount injected in lime kiln")        # Hydrogen - Annual amount injected in lime kiln
    h_to  = m.addVar(name = "Hydrogen - Total amount needed")        # Hydrogen - Total amount needed
    c_s   = m.addVar(name = "CO2 - Annual sold amount")        # CO2 - Annual sold amount
    c_to  = m.addVar(name = "CO2 - Total amount needed")        # CO2 - Total amount needed
    l_s   = m.addVar(name = "Lignin - Annual sold amount")        # Lignin - Annual sold amount
    m_es  = m.addVar(name = "Methane - Annual sold amount")        # Methane - Annual sold amount
    m_eig = m.addVar(name = "Methane - Annual amount injected in gas grid")        # Methane - Annual amount injected in gas grid
    m_eil = m.addVar(name = "Methane - Annual amount injected in lime kiln")        # Methane - Annual amount injected in lime kiln
    m_ls  = m.addVar(name = "Methanol - Annual sold amount")        # Methanol - Annual sold amount
    
    #set objective function
    obj = h_s *PS_HC + h_ig *PS_HI     #revenue hydrogen sale & grid injection
    obj += h_t *D_TKM *FC_T *PP_FT /D_T /100     #savings truck replacement
    obj -= h_t *CC_HT /A_HT     #cost truck replacement
    obj += h_f *D_FKM *FC_F *PP_FF /D_F /100     #savings forklift replacement
    obj -= h_f *CC_HF /A_HF     #cost forklift replacement
    obj -= (h_t *D_TKM *FC_HT /D_T /100 + h_f *D_FKM *FC_HF /D_F /100) *CC_HFI *(1 /A_HFI + CO_HFI)     #cost hydrogen filling station
    obj += h_il * PP_NG     #savings hydrogen injection in lime kiln
    obj -= h_to *CC_EL /PR_HPE /OH_EL /0.08 *(1 /A_EL + CO_EL + CC_SEL /A_SEL)     #cost electrolysis
    obj += o_s * PS_O + o_i *PP_O     #revenue oxygen sale and savings internal oxygen use
    obj += c_s *PS_CO2     #revenue co2 sale
    obj -= c_to * CC_CC * (1 /A_CC + CO_CC)     #cost CCS
    obj += l_s * PS_LI     #revenue lignin sale
    obj -= l_s *CC_LI *(1 /A_LI + CO_LI)     #cost lignoboost
    obj += m_es *PS_MEC + m_eig *PS_MEI + m_eil *PP_NG     #revenue methane sale and grid injection and savings methane inection in lime kiln
    obj -= (m_es + m_eig + m_eil) *CC_ME /1000 /OH_ME *(1 /A_ME + CO_ME)     #cost methanation
    obj += m_ls * PS_ML     #revenue methanol sale
    obj -= m_ls *CC_ML /1000 /OH_ML *(1 /A_ML + CO_ML)     #cost methanol synthesis
    m.setObjective(obj, GRB.MAXIMIZE)
    
    #add constraints
    m.addConstr(h_to == h_s + h_ig + h_t *D_TKM *FC_HT /D_T /100 + h_f *D_FKM *FC_HF /D_F /100 + h_il /39.45 + (m_es + m_eig + m_eil) /PR_MEPH + m_ls /PR_MLPH)
    m.addConstr(c_to == c_s + l_s /PR_LIPCO2 + (m_es + m_eig + m_eil) /PR_MEPCO2 + m_ls /PR_MLPCO2)
    m.addConstr(h_to         <= 1000 *P_EL *PR_HPE *0.08)
    m.addConstr(h_s          <= SP_HC)
    m.addConstr(h_ig         <= SP_HI)
    m.addConstr(h_t          <= D_T)
    m.addConstr(h_f          <= D_F)
    m.addConstr(h_il         <= D_NGLK * R_HLK)
    m.addConstr(o_s + o_i    <= h_to *PR_OPE /PR_HPE /0.08)
    m.addConstr(o_s          <= SP_O)
    m.addConstr(o_i          <= D_O)
    m.addConstr(c_to         <= EM_LK)
    m.addConstr(c_s          <= SP_CO2)
    m.addConstr(l_s          <= SP_LI)
    m.addConstr(m_es         <= SP_MEC)
    m.addConstr(m_eig        <= SP_MEI)
    m.addConstr(h_il + m_eil <= D_NGLK)
    m.addConstr(m_ls         <= SP_ML)
    
    #optimize model
    m.optimize()
    
    for v in m.getVars():
        print('%s %g' % (v.varName, v.x))

    print('Obj: %g' % m.objVal)
        
    return m

def econ_nlin_opt():
    print("-----Economic optimization with non-linear cost structure-----")
    
    #create a new model
    m = gp.Model("co2")
    
    #create decision variables
    o_s   = m.addVar(name = "Oxygen - Annual sold amount")        # Oxygen - Annual sold amount
    o_i   = m.addVar(name = "Oxygen - Annual internally used amount")        # Oxygen - Annual internally used amount
    h_s   = m.addVar(name = "Hydrogen - Annual sold amount")        # Hydrogen - Annual sold amount
    h_ig  = m.addVar(name = "Hydrogen - Annual amount injected in gas grid")        # Hydrogen - Annual amount injected in gas grid
    h_t   = m.addVar(name = "Hydrogen - Number of trucks replaced with fuel cell vehicles", vtype = GRB.INTEGER)        # Hydrogen - Number of trucks replaced with fuel cell vehicles
    h_f   = m.addVar(name = "Hydrogen - Number of forklifts replaced with fuel cell vehicles", vtype = GRB.INTEGER)        # Hydrogen - Number of forklifts replaced with fuel cell vehicles
    h_il  = m.addVar(name = "Hydrogen - Annual amount injected in lime kiln")        # Hydrogen - Annual amount injected in lime kiln
    h_to  = m.addVar(name = "Hydrogen - Total amount needed")        # Hydrogen - Total amount needed
    c_s   = m.addVar(name = "CO2 - Annual sold amount")        # CO2 - Annual sold amount
    c_to  = m.addVar(name = "CO2 - Total amount needed")        # CO2 - Total amount needed
    l_s   = m.addVar(name = "Lignin - Annual sold amount")        # Lignin - Annual sold amount
    m_es  = m.addVar(name = "Methane - Annual sold amount")        # Methane - Annual sold amount
    m_eig = m.addVar(name = "Methane - Annual amount injected in gas grid")        # Methane - Annual amount injected in gas grid
    m_eil = m.addVar(name = "Methane - Annual amount injected in lime kiln")        # Methane - Annual amount injected in lime kiln
    m_ls  = m.addVar(name = "Methanol - Annual sold amount")        # Methanol - Annual sold amount
    
    #gurobi can only solve non-convex problems of quadratic nature
    #other power functions can be piecewise-approximated as linear functions with the Model.addGenConstrPow() function
    #for this, constraints of the format y = x^a need to be introduced
    #to bring the optimization model to this form, some helper variables are defined here
    help_hfi  = m.addVar(name = "help_hfi")
    help_hfi2 = m.addVar(name = "help_hfi2")
    help_el   = m.addVar(name = "help_el")
    help_el2  = m.addVar(name = "help_el2")
    help_cc   = m.addVar(name = "help_cc")
    help_cc2  = m.addVar(name = "help_cc2")
    help_li   = m.addVar(name = "help_li")
    help_li2  = m.addVar(name = "help_li2")
    help_me   = m.addVar(name = "help_me")
    help_me2  = m.addVar(name = "help_me2")
    help_ml   = m.addVar(name = "help_ml")
    help_ml2  = m.addVar(name = "help_ml2")
    
    #set objective function
    obj = h_s *PS_HC + h_ig *PS_HI     #revenue hydrogen sale & grid injection
    obj += h_t *D_TKM *FC_T *PP_FT /D_T /100     #savings truck replacement
    obj -= h_t *CC_HT /A_HT     #cost truck replacement
    obj += h_f *D_FKM *FC_F *PP_FF /D_F /100     #savings forklift replacement
    obj -= h_f *CC_HF /A_HF     #cost forklift replacement 
    obj -= CC_KHFI * help_hfi *(1 /A_HFI + CO_HFI)     #cost hydrogen filling station
    obj += h_il * PP_NG     #savings hydrogen injection in lime kiln
    obj -= CC_KEL *help_el *(1 /A_EL + CO_EL + CC_SEL /A_SEL)     #cost electrolysis
    obj += o_s * PS_O + o_i *PP_O     #revenue oxygen sale and savings internal oxygen use
    obj += c_s *PS_CO2     #revenue co2 sale
    obj -= CC_KCC *help_cc * (1 /A_CC + CO_CC)     #cost CCS
    obj += l_s * PS_LI     #revenue lignin sale
    obj -= CC_KLI *help_li *(1 /A_LI + CO_LI)     #cost lignoboost
    obj += m_es *PS_MEC + m_eig *PS_MEI + m_eil *PP_NG     #revenue methane sale and grid injection and savings methane inection in lime kiln    
    obj -= CC_KME *help_me *(1 /A_ME + CO_ME)     #cost methanation
    obj += m_ls * PS_ML     #revenue methanol sale
    obj -= CC_KML *help_ml *(1 /A_ML + CO_ML)     #cost methanol synthesis
    m.setObjective(obj, GRB.MAXIMIZE)
    
    #add constraints
    m.addConstr(h_to == h_s + h_ig + h_t *D_TKM *FC_HT /D_T /100 + h_f *D_FKM *FC_HF /D_F /100 + h_il /39.45 + (m_es + m_eig + m_eil) /PR_MEPH + m_ls /PR_MLPH)
    m.addConstr(c_to == c_s + l_s /PR_LIPCO2 + (m_es + m_eig + m_eil) /PR_MEPCO2 + m_ls /PR_MLPCO2)
    m.addConstr(h_to         <= 1000 *P_EL *PR_HPE *0.08)
    m.addConstr(h_s          <= SP_HC)
    m.addConstr(h_ig         <= SP_HI)
    m.addConstr(h_t          <= D_T)
    m.addConstr(h_f          <= D_F)
    m.addConstr(h_il         <= D_NGLK * R_HLK)
    m.addConstr(o_s + o_i    <= h_to *PR_OPE /PR_HPE /0.08)
    m.addConstr(o_s          <= SP_O)
    m.addConstr(o_i          <= D_O)
    m.addConstr(c_to         <= EM_LK)
    m.addConstr(c_s          <= SP_CO2)
    m.addConstr(l_s          <= SP_LI)
    m.addConstr(m_es         <= SP_MEC)
    m.addConstr(m_eig        <= SP_MEI)
    m.addConstr(h_il + m_eil <= D_NGLK)
    m.addConstr(m_ls         <= SP_ML)
    
    #add helper constraints for linear approximations of non-convex parts of the objective function
    m.addConstr(help_hfi2 == ((h_t *D_TKM *FC_HT /D_T /100 + h_f *D_FKM *FC_HF /D_F /100) /CA_KHFI))
    m.addGenConstrPow(help_hfi2, help_hfi, N_KHFI)
    m.addConstr(help_el2 == (h_to /CA_KEL /PR_HPE /OH_EL /0.08))
    m.addGenConstrPow(help_el2, help_el, N_KEL)
    m.addConstr(help_cc2 == (c_to /CA_KCC))
    m.addGenConstrPow(help_cc2, help_cc, N_KCC)
    m.addConstr(help_li2 == (l_s /CA_KLI))
    m.addGenConstrPow(help_li2, help_li, N_KLI)
    m.addConstr(help_me2 == ((m_es + m_eig + m_eil) /CA_KME /1000 /OH_ME))
    m.addGenConstrPow(help_me2 ,help_me ,N_KME)
    m.addConstr(help_ml2 == (m_ls /CA_KML /1000 /OH_ML))
    m.addGenConstrPow(help_ml2, help_ml, N_KML)
    
    #optimize model
    m.optimize()
    
    for v in m.getVars():
        print('%s %g' % (v.varName, v.x))

    print('Obj: %g' % m.objVal)
        
    return m

#calculates the CO2 result based on a given optimzation result
def co2_fct(opt):
    
    #read variables from optimzation result
    var = opt.getVars()
    o_s   = var[0].x        # Oxygen - Annual sold amount
    o_i   = var[1].x        # Oxygen - Annual internally used amount
    h_s   = var[2].x        # Hydrogen - Annual sold amount
    h_ig  = var[3].x        # Hydrogen - Annual amount injected in gas grid
    h_t   = var[4].x        # Hydrogen - Number of trucks replaced with fuel cell vehicles
    h_f   = var[5].x        # Hydrogen - Number of forklifts replaced with fuel cell vehicles
    h_il  = var[6].x        # Hydrogen - Annual amount injected in lime kiln
    h_to  = var[7].x        # Hydrogen - Total amount needed
    c_s   = var[8].x        # CO2 - Annual sold amount
    c_to  = var[9].x        # CO2 - Total amount needed
    l_s   = var[10].x       # Lignin - Annual sold amount
    m_es  = var[11].x       # Methane - Annual sold amount
    m_eig = var[12].x       # Methane - Annual amount injected in gas grid
    m_eil = var[13].x       # Methane - Annual amount injected in lime kiln
    m_ls  = var[14].x       # Methanol - Annual sold amount
    
    #calculate objective function
    #obj = EM_T * h_t /1000000
    #obj += EM_F * h_f /1000000
    obj =  h_t * EM_T * D_TKM / 1000000 / D_T
    obj += h_f * EM_F * D_FKM / 1000000 / D_F
    obj += EM_LK / D_NGLK * h_il
    obj += c_s
    
    return obj
    
#calculates the linear economic result based on a given optimzation result
def econ_lin_fct(opt):
    
    #read variables from optimzation result
    var = opt.getVars()
    o_s   = var[0].x        # Oxygen - Annual sold amount
    o_i   = var[1].x        # Oxygen - Annual internally used amount
    h_s   = var[2].x        # Hydrogen - Annual sold amount
    h_ig  = var[3].x        # Hydrogen - Annual amount injected in gas grid
    h_t   = var[4].x        # Hydrogen - Number of trucks replaced with fuel cell vehicles
    h_f   = var[5].x        # Hydrogen - Number of forklifts replaced with fuel cell vehicles
    h_il  = var[6].x        # Hydrogen - Annual amount injected in lime kiln
    h_to  = var[7].x        # Hydrogen - Total amount needed
    c_s   = var[8].x        # CO2 - Annual sold amount
    c_to  = var[9].x        # CO2 - Total amount needed
    l_s   = var[10].x       # Lignin - Annual sold amount
    m_es  = var[11].x       # Methane - Annual sold amount
    m_eig = var[12].x       # Methane - Annual amount injected in gas grid
    m_eil = var[13].x       # Methane - Annual amount injected in lime kiln
    m_ls  = var[14].x       # Methanol - Annual sold amount
    
    #calculate objective function
    obj = h_s *PS_HC + h_ig *PS_HI     #revenue hydrogen sale & grid injection
    obj += h_t *D_TKM *FC_T *PP_FT /D_T /100     #savings truck replacement
    obj -= h_t *CC_HT /A_HT     #cost truck replacement
    obj += h_f *D_FKM *FC_F *PP_FF /D_F /100     #savings forklift replacement
    obj -= h_f *CC_HF /A_HF     #cost forklift replacement
    obj -= (h_t *D_TKM *FC_HT /D_T /100 + h_f *D_FKM *FC_HF /D_F /100) *CC_HFI *(1 /A_HFI + CO_HFI)     #cost hydrogen filling station
    obj += h_il * PP_NG     #savings hydrogen injection in lime kiln
    obj -= h_to *CC_EL /PR_HPE /OH_EL /0.08 *(1 /A_EL + CO_EL + CC_SEL /A_SEL)     #cost electrolysis
    obj += o_s * PS_O + o_i *PP_O     #revenue oxygen sale and savings internal oxygen use
    obj += c_s *PS_CO2     #revenue co2 sale
    obj -= c_to * CC_CC * (1 /A_CC + CO_CC)     #cost CCS
    obj += l_s * PS_LI     #revenue lignin sale
    obj -= l_s *CC_LI *(1 /A_LI + CO_LI)     #cost lignoboost
    obj += m_es *PS_MEC + m_eig *PS_MEI + m_eil *PP_NG     #revenue methane sale and grid injection and savings methane inection in lime kiln
    obj -= (m_es + m_eig + m_eil) *CC_ME /1000 /OH_ME *(1 /A_ME + CO_ME)     #cost methanation
    obj += m_ls * PS_ML     #revenue methanol sale
    obj -= m_ls *CC_ML /1000 /OH_ML *(1 /A_ML + CO_ML)     #cost methanol synthesis
    
    return obj

#calculates the non-linear economic result based on a given optimzation result
def econ_nlin_fct(opt):
    
    #read variables from optimzation result
    var = opt.getVars()
    o_s   = var[0].x        # Oxygen - Annual sold amount
    o_i   = var[1].x        # Oxygen - Annual internally used amount
    h_s   = var[2].x        # Hydrogen - Annual sold amount
    h_ig  = var[3].x        # Hydrogen - Annual amount injected in gas grid
    h_t   = var[4].x        # Hydrogen - Number of trucks replaced with fuel cell vehicles
    h_f   = var[5].x        # Hydrogen - Number of forklifts replaced with fuel cell vehicles
    h_il  = var[6].x        # Hydrogen - Annual amount injected in lime kiln
    h_to  = var[7].x        # Hydrogen - Total amount needed
    c_s   = var[8].x        # CO2 - Annual sold amount
    c_to  = var[9].x        # CO2 - Total amount needed
    l_s   = var[10].x       # Lignin - Annual sold amount
    m_es  = var[11].x       # Methane - Annual sold amount
    m_eig = var[12].x       # Methane - Annual amount injected in gas grid
    m_eil = var[13].x       # Methane - Annual amount injected in lime kiln
    m_ls  = var[14].x       # Methanol - Annual sold amount
    
    #calculate helper variables from optimzation result
    help_hfi2 = ((h_t *D_TKM *FC_HT /D_T /100 + h_f *D_FKM *FC_HF /D_F /100) /CA_KHFI)
    help_hfi  = help_hfi2**N_KHFI
    help_el2  = (h_to /CA_KEL /PR_HPE /OH_EL /0.08)
    help_el   = help_el2**N_KEL
    help_cc2  = (c_to /CA_KCC)
    help_cc   = help_cc2**N_KCC
    help_li2  = (l_s /CA_KLI)
    help_li   = (help_li2**N_KLI)
    help_me2  = ((m_es + m_eig + m_eil) /CA_KME /1000 /OH_ME)
    help_me   = help_me2**N_KME
    help_ml2  = (m_ls /CA_KML /1000 /OH_ML)
    help_ml   = help_ml2**N_KML

    #calculate objective function
    obj = h_s *PS_HC + h_ig *PS_HI     #revenue hydrogen sale & grid injection
    obj += h_t *D_TKM *FC_T *PP_FT /D_T /100     #savings truck replacement
    obj -= h_t *CC_HT /A_HT     #cost truck replacement
    obj += h_f *D_FKM *FC_F *PP_FF /D_F /100     #savings forklift replacement
    obj -= h_f *CC_HF /A_HF     #cost forklift replacement 
    obj -= CC_KHFI * help_hfi *(1 /A_HFI + CO_HFI)     #cost hydrogen filling station
    obj += h_il * PP_NG     #savings hydrogen injection in lime kiln
    obj -= CC_KEL *help_el *(1 /A_EL + CO_EL + CC_SEL /A_SEL)     #cost electrolysis
    obj += o_s * PS_O + o_i *PP_O     #revenue oxygen sale and savings internal oxygen use
    obj += c_s *PS_CO2     #revenue co2 sale
    obj -= CC_KCC *help_cc * (1 /A_CC + CO_CC)     #cost CCS
    obj += l_s * PS_LI     #revenue lignin sale
    obj -= CC_KLI *help_li *(1 /A_LI + CO_LI)     #cost lignoboost
    obj += m_es *PS_MEC + m_eig *PS_MEI + m_eil *PP_NG     #revenue methane sale and grid injection and savings methane inection in lime kiln    
    obj -= CC_KME *help_me *(1 /A_ME + CO_ME)     #cost methanation
    obj += m_ls * PS_ML     #revenue methanol sale
    obj -= CC_KML *help_ml *(1 /A_ML + CO_ML)     #cost methanol synthesis
    
    return obj

#---------------------------------------------------------------------------

#call of optimization
opt(co2=co2, econ_lin=econ_lin, econ_nlin=econ_nlin)
#econ_nlin_fct(opt)