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.Obj(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]))