from gekko import GEKKO

#  Initialize Gekko model
m = GEKKO()

# Optimal heat pump and regenerative exchanger that minimizes the total present
#  worth of costs (capital cost and operating cost). Specifically, determine the
#  areas of the heat exchangers (evaporator, condensers, and regenerative
#  exchanger), size of the compressor, and temperatures t1, t2, te and tc 
#  that result in the minimum total cost

# Parameters
U = 0.6
Ureg = 0.5
mdotm = 4
cpm = 3.75
cpw = 1.00
intrate = 0.09
nper = 6

# Inlet temperature (degC)
ti = 7              

# Outlet temperature (degC)
to = 4

# Pasteurization Temperature (degC), target temperature
tp = 73

# Variables
# compressor work (kW)
W = m.Var()

# heat exchanger areas (m^2)
Ae = m.Var(value = 10, lb = 1.0, ub = 1000)
Areg = m.Var(value = 10, lb = 1.0, ub = 1000)
Afc = m.Var(value = 10, lb = 1.0, ub = 1000)
Aac = m.Var(value = 10, lb = 1.0, ub = 1000)

t1 = m.Var(value = ti)             # Temperature 1 (degC)
t2 = m.Var(value = ti)             # Temperature 2 (degC)
te = m.Var(value = ti)             # Evaporator Temperature (degC)
tc = m.Var(value = ti)             # Condenser Temperature (degC)
Ctot = m.Var(value = 100000)       # Total Cost

# evaporator
Qe  = m.Var(value = 1)
dapp1e = m.Var(value = 1)
dapp2e = m.Var(value = 1)

# regenerative heat exchanger 
Qreg = m.Var(value = 1)
dapp1reg = m.Var(value = 1)
dapp2reg = m.Var(value = 1)

# fore condenser 
Qfc = m.Var(value = 1)
dapp1fc = m.Var(value = 1)
dapp2fc = m.Var(value = 1)

# after condenser 
Qac = m.Var(value = 1)
dapp1ac = m.Var(value = 1)
dapp2ac = m.Var(value = 1)

# temps of approach 
dappe = m.Var(lb = 0.1)     # >= 10
dappc = m.Var(lb = 0.1)     # >= 10
dappreg1 = m.Var(lb = 0.1)  # >= 10
dappreg2 = m.Var(lb = 0.1)  # >= 10

# ratios
ratio_e = m.Var(value = .05, lb = 0.001, ub = 0.999)
ratio_reg = m.Var(value = .05, lb = 0.000, ub = 1.0)
ratio_fc = m.Var(value = .05, lb = 0.001, ub = 0.999)
ratio_ac = m.Var(value = .05, lb = 0.001, ub = 0.999)

# Log-mean temperature differences
delte = m.Var(value = 1)
deltreg = m.Var(value = 1)
deltfc = m.Var(value = 1)
deltac = m.Var(value = 1)

# Coefficient Of Performance
COP = m.Var(value = 1)

# Intermediates
# Cost elec year at 4 hours per day
Celyear= m.Intermediate(W * 4.0 * 365.0 * 0.07)

# Cost total over time
# Multiplier is about 4.5 (instead of 6.0) because of the 
# interest rate on unspent capital
mult = m.Intermediate((((1.0+intrate)**nper)-1.0) / (intrate*(1.0+intrate)**nper))
Celtot = m.Intermediate(Celyear * mult)

# Cost equipment
Cequip = m.Intermediate(Ae*200.0 +  Afc*200.0 + Aac*200.0 + Areg*200 + W*240.0)

#temperatures in Kelvin
TeK = m.Intermediate(te + 273)
TcK = m.Intermediate(tc + 273)

# Equations
m.Equations([
    # Total cost
    Ctot == Celtot + Cequip,

    # overall energy balance on regen heat exchanger
    # mdotm and cpm cancel from each term
    t1 + t2 == tp + ti,

    # log-mean temperature difference
    #   lmtd(thi,tco,tho,tci)

    # evaporator 
    Qe == mdotm * cpm *(t2 - to),
    dapp1e == to - te,
    dapp2e == t2 - te,

    # regenerative heat exchanger 
    Qreg == mdotm * cpm * (tp-t2),
    dapp1reg == t2 - ti,
    dapp2reg == tp - t1,

    # fore condenser 
    Qfc == mdotm * cpm * (tp - t1),
    dapp1fc == tc - tp,
    dapp2fc == tc - t1,

    # after condenser 
    Qac == W + mdotm * cpm * (ti - to),
    dapp1ac == tc - 35.0,
    dapp2ac == tc - 30.0,

    # approach temperatures 
    dappe == to - te,         # Temp app evap
    dappc == tc - tp,         # Temp app fore cond
    dappreg1 == t2 - ti,      # Temp app reg1
    dappreg2 == tp - t1,      # Temp app reg2

    # ratios
    ratio_e   * dapp2e   == dapp1e, 
    ratio_reg * dapp2reg == dapp1reg,
    ratio_fc  * dapp2fc  == dapp1fc,
    ratio_ac  * dapp2ac  == dapp1ac,

    # Log-mean Temperature Difference
    delte   * m.log(ratio_e)   == (dapp1e - dapp2e),
    #deltreg * m.log(ratio_reg) = (dapp1reg - dapp2reg)
    deltfc  * m.log(ratio_fc)  == (dapp1fc - dapp2fc),
    deltac  * m.log(ratio_ac)  == (dapp1ac - dapp2ac),

    # Average temperature difference
    #delte    = (dapp1e   + dapp2e)   / 2
    deltreg  == (dapp1reg + dapp2reg) / 2,
    #deltfc   = (dapp1fc  + dapp2fc)  / 2
    #deltac   = (dapp1ac  + dapp2ac)  / 2

    # coefficient of performance
    # 75% of Carnot Efficiency
    COP * (TcK-TeK) == 0.75 * TeK,

    # Work to compressor 
    W * COP == Qe,

    # heat transferred for each exchanger
    Ae   * U    * delte   == Qe,
    Areg * Ureg * deltreg == Qreg,
    Afc  * U    * deltfc  == Qfc,
    Aac  * U    * deltac  == Qac
    ])

# objective function
m.Minimize(Ctot)

# minimize ctot
m.options.SOLVER = 1
m.options.IMODE = 3

m.solve()

print('Minimum cost: ' + str(Ctot[0]))