from gekko import GEKKO

m = GEKKO()

#%% Model Oxygen Supply
#Constants
R = 670.6     #gas law (cubic feet-psia/(HP - HR))
T = 530       #Temperature  (degrees Rankine)
z = 0.98      #Compressibility factor ()
M = 15.999    #Molecular weight (u)
k_1 = 14005.8 #unit conversion factor

#cycle parameters
N = 8000    #Number of cycles per year
D_0 = 2.4   #Low demand rate of O_2 (tons/hr)
D_1 = 37    #High demand rate of O_2 (tons/hr)
t_1 = 0.6   #time interval of high demand rate
t_2 = 1     #time interval of one cycle of low and high demand rate

#cost parameters
a_1 = 60.9
a_2 = 5.83
b_1 = 2.5e-5
b_2 = 680      
b_3 = 0.85
b_4 = 6.0e-3   #cost of power
s_1 = 3.0e-5
s_2 = 374
s_3 = 0.9
d = 5           #annual cost factor

#physical parameters
p_0 = 200       #O_2 delivery pressure
k_2 = 0.75      #Compressor efficiency

#%% Variables
#Independent variables
F = m.Var()             #Production rate(tons/hr)
p = m.Var()             #Tank pressure (psia)
annual_cost = m.Var()   #Objective. 

#%% Intermediates
I_max = m.Intermediate((D_1 - F)*(t_2 - t_1))
V = m.Intermediate((I_max/M)*(R*T/p)*z)
H = m.Intermediate((I_max/t_1)*(R*T/(k_1*k_2))*m.log(p/p_0))
#Oxygen plant annual cost
C_1 = m.Intermediate(a_1 + a_2*F)
#Capital cost of storage vessels
C_2 = m.Intermediate(s_1*s_2*V**s_3)
#Capital cost of compressors
C_3 = m.Intermediate(0.2*b_1*b_2*H**b_3)
#Compressor power cost
C_4 = m.Intermediate(0.001*b_4*t_1*H)

#%% Equations
m.Equations([
        annual_cost == C_1 + d*(C_2 + C_3) + N*C_4,
        F >= (D_0*t_1 + D_1*(t_2 - t_1))/t_2,
        p >= p_0 + 1
        ])

#%% Objective
m.Minimize(annual_cost)

#%% Solve
m.solve()

print("Optimal annual cost:" , str(annual_cost[0]))
print("Optimal production rate:" , str(F[0]))
print("Optimal maximum tank pressure:" , str(p[0]))