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]))