I am currently in the process of trying to determine the optimal inflow conditions for the refueling process of a hydrogen (H2) vehicle using GEKKO. Below are the coupled, ordinary differential equations which govern how the temperature of H2 and the fuel tank wall change over the fueling time.
T.dt() = (1+alpha)*(T_star - T)/(t_star + t)
T_w.dt() = (T - T_w)/t_w_star
where
alpha = (a_in*A_in)/(c_v*m_dot_in), t_star = m_0/m_dot_in, t_w_star = (m_w*c_w)/(a_in*A_in)
T_star = gamma_p*T_inf + alpha_p*T_w, gamma_p = gamma/(1 + alpha), alpha_p = alpha/(1 + alpha)
Here, m_0 is the initial mass of H2 in the tank, m_dot_in is the mass flow rate of H2 into the tank, gamma is the ratio of the specific heats for H2, T_inf is the inflow temperature of H2, and the other variables are intermediate variables/tank parameters. Through the refueling process, m_dot_in is taken to be constant (but unknown), so the mass of H2 in the tank over time is defined as:
m = m_0 + m_dot_in*t
Additionally, the pressure of H2 within the tank may be calculated with a real gas equation of state (I use the Peng-Robinson equation of state for this model).
What I am trying to do with this model is determine the optimal m_dot_in, T_inf, and m_0 to minimize the total fueling time, t_f. Some constraints on the variables are that T<=358.15 K throughout the whole refueling process (for safety reasons), and that the final pressure of H2 within the tank must be 35 MPa. For this model, I consider t_f, m_dot_in, m_0, and T_inf to be fixed variables that have the following bounds:
60 sec <= t_f <= 300 sec
0.0005 kg/sec <= m_dot_in <= 0.03 kg/sec
5% of m_f <= m_0 <= 90% of m_f, where m_f = 1.79 kg
288.15 K <= T_inf <= 303.15 K
Below I have copied my code for this optimization problem using GEKKO:
# HYDROGEN TANK REFUELING MODEL
# OPTIMIZE MODEL WITH GEKKO OPTIMIZATION SUITE
from gekko import GEKKO
m = GEKKO()
# CONSTANTS
## TANK PARAMETERS (ASSUME TYPE III, ALUMINUM, 74 l, RATED FOR 35 MPa)
V = 0.074 # m^3
a_in = 167/19e-3 # W/m^2/K
c_w = 2730 # J/kg/K
rho_w = 900 # kg/m^3
m_w = rho_w*(np.pi*(((0.358+19e-3)/2)**2)*(0.893+19e-3) - np.pi*((0.358/2)**2)*0.893)
A_in = 2*np.pi*(0.358/2)*((0.358/2) + 0.893) # m^2
T_w0 = 293.15 # K
m_f = 1.79 # final mass of hydrogen in tank, kg
## HYDROGEN PARAMETERS
c_p = 14.615e3 # specific heat at constant pressure, J/kg/K
c_v = 10.316e3 # specific heat at constant volume, J/kg/K
gamma = c_p/c_v
R = 8.314/M_H2 # gas constant for H2, J/kgK
T_c = -240 + 273.15 # critical temperature for H2, K
p_c = 1.3e6 # critical pressure for H2, Pa
w_H2 = -0.219 # acentric factor for H2
a = 0.45724*(R**2 * T_c**2)/(p_c**2)
b = 0.0778*(R*T_c)/p_c
kappa = 0.37464 + 1.54226*w_H2 - 0.26992*(w_H2**2)
# SET TIME ANALYSIS POINTS
nt = 101
tm = np.linspace(0, 1, nt)
m.time = tm
# options
m.options.NODES = 6
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.MAX_ITER = 500
m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0
# FIXED VARIABLES
t_f = m.FV(value=60.0,lb=60.0,ub=300.0) # final fuel time, s
t_f.STATUS = 1
m_dot_in = m.FV(value=0.001,lb=0.0005,ub=0.03) # mass flow rate into tank, kg/s
m_dot_in.STATUS = 1
m_0 = m.FV(value=0.1*m_f,lb=0.05*m_f,ub=0.9*m_f) # initial mass of H2 in tank (as % of m_f), kg
m_0.STATUS = 1
T_inf = m.FV(value=20 + 273.15,lb=15 + 273.15,ub=30 + 273.15) # inflow temperature, K
T_inf.STATUS = 1
# PARAMETERS
f = np.zeros(nt)
f[-1] = 1.0
final = m.Param(value=f)
# VARIABLES
T = m.Var(value=15+273.15,lb=15+273.15,ub=85+273.15)
T_w = m.Var(value=T_w0,lb=T_w0,ub=85+273.15)
mass = m.Var(value=m_0,lb=m_0,ub=m_f)
p = m.Var(value=1.0e6,lb=0.0,ub=35.0e6)
# INTERMEDIATES
alpha = m.Intermediate((a_in*A_in)/c_v/m_dot_in)
gamma_p = m.Intermediate(gamma/(1 + alpha))
alpha_p = m.Intermediate(alpha/(1 + alpha))
t_star = m.Intermediate(m_0/m_dot_in)
t_w_star = m.Intermediate((m_w*c_w)/(a_in*A_in))
T_star = m.Intermediate(gamma_p*T_inf + alpha_p*T_w)
alpha_T = m.Intermediate(1 + kappa*(1 - (T/T_c)**0.5))
v = m.Intermediate(V/mass) # specific volume, m^3/kg
# EQUATIONS
m.Equation(mass==t_f*(m_0 + m_dot_in*m.time))
m.Equation(T.dt()==t_f*(1 + alpha)*((T_star-T)/(t_star+tm)))
m.Equation(T_w.dt()==t_f*((T-T_w)/t_w_star))
m.Equation(p*1.0e6==t_f*((R*T/(v-b)) - ((a*alpha_T**2)/(v*(v+b) + b*(v-b)))))
m.Equation(T*final<=85+273.15)
m.Equation(T_w*final<=85+273.15)
# SPECIFIY ENDPOINT CONDITIONS
m.fix(mass, pos=len(m.time)-1, val=m_f)
m.fix(p, pos=len(m.time)-1, val=35.0e6)
# MINIMIZE FINAL FUEL TIME
m.Obj(t_f)
# SOLVE
m.solve()
# RESULTS
print('Final Time: ' + str(t_f.value[0]))
This code currently gives me the following error:
apm 45.3.69.90_gk_model46 <br><pre> ----------------------------------------------------------------
APMonitor, Version 1.0.1
APMonitor Optimization Suite
----------------------------------------------------------------
@error: Equation Definition
Equation without an equality (=) or inequality (>,<)
0.140.150.160.170.180.190.20.210.220.230.240.250.260.27
STOPPING...
---------------------------------------------------------------------------
Exception Traceback (most recent call last)
<ipython-input-102-4d40bf2f7c9c> in <module>
87
88 # SOLVE
---> 89 m.solve()
90
91 # RESULTS
~\anaconda3\lib\site-packages\gekko\gekko.py in solve(self, disp, debug, GUI, **kwargs)
2172 #print APM error message and die
2173 if (debug >= 1) and ('@error' in response):
-> 2174 raise Exception(response)
2175
2176 #load results
Exception: @error: Equation Definition
Equation without an equality (=) or inequality (>,<)
0.140.150.160.170.180.190.20.210.220.230.240.250.260.27
STOPPING...
I am very new to optimization in general, and I tried including several different equality and inequality constraints, but nothing seems to work. I thought I was doing it correctly based off example problems and information from the APMonitor website, but obviously something is off with my implementation. I was wondering if someone knew what I should change/add or if I am doing something completely wrong? Any help would be much appreciated!
Thank you for your time,
Evan
EDIT: Based on Dr. Hedengren's answer, I tried to simplifying the model such that the variables of mass and p were not included, as they only depend on the final values for t_f, m_dot_in, and T and may be calculated after the solution has been obtained. Below is my edited code:
# HYDROGEN TANK REFUELING MODEL
# OPTIMIZE MODEL WITH GEKKO OPTIMIZATION SUITE
from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
# CONSTANTS
## TANK PARAMETERS (ASSUME TYPE III, ALUMINUM, 74 l, RATED FOR 35 MPa)
V = 0.074 # m^3
a_in = 167/19e-3 # W/m^2/K
c_w = 2730 # J/kg/K
rho_w = 900 # kg/m^3
m_w = rho_w*(np.pi*(((0.358+19e-3)/2)**2)*(0.893+19e-3) - np.pi*((0.358/2)**2)*0.893)
A_in = 2*np.pi*(0.358/2)*((0.358/2) + 0.893) # m^2
T_w0 = 293.15 # K
m_f = 1.79 # final mass of hydrogen in tank, kg
## HYDROGEN PARAMETERS
c_p = 14.615e3 # specific heat at constant pressure, J/kg/K
c_v = 10.316e3 # specific heat at constant volume, J/kg/K
gamma = c_p/c_v
M_H2 = 2.02
R = 8.314/M_H2 # gas constant for H2, J/kgK
T_c = -240 + 273.15 # critical temperature for H2, K
p_c = 1.3e6 # critical pressure for H2, Pa
w_H2 = -0.219 # acentric factor for H2
a = 0.45724*(R**2 * T_c**2)/(p_c**2)
b = 0.0778*(R*T_c)/p_c
kappa = 0.37464 + 1.54226*w_H2 - 0.26992*(w_H2**2)
## SET INFLOW TEMPERATURE AND INITIAL MASS IN TANK
m_0 = 0.1*m_f
T_inf = 20 + 273.15
# SET TIME ANALYSIS POINTS
nt = 101
tm = np.linspace(0, 1, nt)
m.time = tm
t = m.Param(tm, name='time')
# options
m.options.NODES = 6
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.MAX_ITER = 500
m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0
# FIXED VARIABLES
t_f = m.FV(value=60.0,lb=60.0,ub=300.0) # final fuel time, s
t_f.STATUS = 0
m_dot_in = m.FV(value=0.001,lb=0.0005,ub=0.03) # mass flow rate into tank, kg/s
m_dot_in.STATUS = 0
# m_0 = m.FV(value=0.1*m_f,lb=0.05*m_f,ub=0.9*m_f) # initial mass of H2 in tank (as % of m_f), kg
# m_0.STATUS = 0
# T_inf = m.FV(value=20 + 273.15,lb=15 + 273.15,ub=30 + 273.15) # inflow temperature, K
# T_inf.STATUS = 0
# PARAMETERS
f = np.zeros(nt)
f[-1] = 1.0
final = m.Param(value=f, name='final')
# VARIABLES
T = m.Var(value=15+273.15,lb=15+273.15,ub=85+273.15, name='H2 Temp')
T_w = m.Var(value=T_w0,lb=T_w0,ub=85+273.15, name='Wall Temp')
# mass = m.Var(value=m_0,lb=m_0,ub=m_f, name='H2 Mass')
# p = m.Var(value=1.0e6,lb=0.0,ub=35.0e6, name='H2 Press')
# INTERMEDIATES
alpha = m.Intermediate((a_in*A_in)/c_v/m_dot_in, name='alpha')
gamma_p = m.Intermediate(gamma/(1 + alpha), name='gamma_p')
alpha_p = m.Intermediate(alpha/(1 + alpha), name='alpha_p')
t_star = m.Intermediate(m_0/m_dot_in, name='t_star')
t_w_star = m.Intermediate((m_w*c_w)/(a_in*A_in), name='t_w_star')
T_star = m.Intermediate(gamma_p*T_inf + alpha_p*T_w, name='Temp_star')
# alpha_T = m.Intermediate(1 + kappa*(1 - (T/T_c)**0.5))
# v = m.Intermediate(V/mass) # specific volume, m^3/kg
# EQUATIONS
# m.Equation(mass==t_f*(m_0 + m_dot_in*t*t_f))
m.Equation(T.dt()==t_f*(1 + alpha)*((T_star-T)/(t_star+t*t_f)))
m.Equation(T_w.dt()==t_f*((T-T_w)/t_w_star))
# m.Equation(p==t_f*((R*T/(v-b)) - ((a*alpha_T**2)/(v*(v+b) + b*(v-b)))))
# m.Equation((T-(85+273.15))*final<=0)
# m.Equation((T_w-(85+273.15))*final<=0)
# SPECIFIY ENDPOINT CONDITIONS
# m.Minimize(final*(mass-m_f)**2)
# m.Minimize(final*(p-35.0e6)**2)
m.Minimize(final*(T-351)**2)
#m.fix(mass, pos=len(m.time)-1, val=m_f)
#m.fix(p, pos=len(m.time)-1, val=35.0e6)
# MINIMIZE FINAL FUEL TIME
m.Minimize(t_f)
# SOLVE
m.options.SOLVER = 3
m.open_folder()
m.solve()
# RESULTS
print('Final Time: ' + str(t_f.value[0]))
I am still getting infeasibilities (not as many as before), but I am having trouble understanding what said infeasibilities mean and how to go about fixing them. Below are the infeasibilities I am getting:
************************************************
***** POSSIBLE INFEASBILE EQUATIONS ************
************************************************
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
1 0.0000E+00 -7.5600E-04 0.0000E+00 7.5600E-04 p(1).n(2).Eqn(1): 0 = $h2_temp-(((((p2)*((1+alpha))))*((((temp_star-h2_temp))/((t_star+((time)*(p2))))))))
Variable Lower Value Upper $Value Name
0 -1.2346E+20 1.0000E+00 1.2346E+20 0.0000E+00 p(1).n(2).time
0 6.0000E+01 6.0000E+01 3.0000E+02 0.0000E+00 p(1).n(1).p2
0 5.0000E-04 1.0000E-03 3.0000E-02 0.0000E+00 p(1).n(1).p3
1 2.8815E+02 2.8949E+02 3.5815E+02 9.7624E+02 p(1).n(2).h2_temp
2 2.9315E+02 2.9315E+02 3.5815E+02 -7.9554E+01 p(1).n(2).wall_temp
1 2.8815E+02 2.8949E+02 3.5815E+02 9.7624E+02 p(1).n(2).h2_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
5 0.0000E+00 -7.5600E-04 0.0000E+00 7.5600E-04 p(1).n(3).Eqn(1): 0 = $h2_temp-(((((p2)*((1+alpha))))*((((temp_star-h2_temp))/((t_star+((time)*(p2))))))))
Variable Lower Value Upper $Value Name
0 -1.2346E+20 1.0000E+00 1.2346E+20 0.0000E+00 p(1).n(3).time
0 6.0000E+01 6.0000E+01 3.0000E+02 0.0000E+00 p(1).n(1).p2
0 5.0000E-04 1.0000E-03 3.0000E-02 0.0000E+00 p(1).n(1).p3
3 2.8815E+02 2.9123E+02 3.5815E+02 5.2535E+02 p(1).n(3).h2_temp
4 2.9315E+02 2.9315E+02 3.5815E+02 -4.1620E+01 p(1).n(3).wall_temp
3 2.8815E+02 2.9123E+02 3.5815E+02 5.2535E+02 p(1).n(3).h2_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
9 0.0000E+00 -7.5600E-04 0.0000E+00 7.5600E-04 p(1).n(4).Eqn(1): 0 = $h2_temp-(((((p2)*((1+alpha))))*((((temp_star-h2_temp))/((t_star+((time)*(p2))))))))
Variable Lower Value Upper $Value Name
0 -1.2346E+20 1.0000E+00 1.2346E+20 0.0000E+00 p(1).n(4).time
0 6.0000E+01 6.0000E+01 3.0000E+02 0.0000E+00 p(1).n(1).p2
0 5.0000E-04 1.0000E-03 3.0000E-02 0.0000E+00 p(1).n(1).p3
5 2.8815E+02 2.9229E+02 3.5815E+02 2.5164E+02 p(1).n(4).h2_temp
6 2.9315E+02 2.9315E+02 3.5815E+02 -1.8591E+01 p(1).n(4).wall_temp
5 2.8815E+02 2.9229E+02 3.5815E+02 2.5164E+02 p(1).n(4).h2_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
13 0.0000E+00 -7.5600E-04 0.0000E+00 7.5600E-04 p(1).n(5).Eqn(1): 0 = $h2_temp-(((((p2)*((1+alpha))))*((((temp_star-h2_temp))/((t_star+((time)*(p2))))))))
Variable Lower Value Upper $Value Name
0 -1.2346E+20 1.0000E+00 1.2346E+20 0.0000E+00 p(1).n(5).time
0 6.0000E+01 6.0000E+01 3.0000E+02 0.0000E+00 p(1).n(1).p2
0 5.0000E-04 1.0000E-03 3.0000E-02 0.0000E+00 p(1).n(1).p3
7 2.8815E+02 2.9274E+02 3.5815E+02 1.3550E+02 p(1).n(5).h2_temp
8 2.9315E+02 2.9315E+02 3.5815E+02 -8.8200E+00 p(1).n(5).wall_temp
7 2.8815E+02 2.9274E+02 3.5815E+02 1.3550E+02 p(1).n(5).h2_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
25 0.0000E+00 -7.5600E-04 0.0000E+00 7.5600E-04 p(2).n(3).Eqn(1): 0 = $h2_temp-(((((p2)*((1+alpha))))*((((temp_star-h2_temp))/((t_star+((time)*(p2))))))))
Variable Lower Value Upper $Value Name
0 -1.2346E+20 1.0000E+00 1.2346E+20 0.0000E+00 p(2).n(3).time
0 6.0000E+01 6.0000E+01 3.0000E+02 0.0000E+00 p(1).n(1).p2
0 5.0000E-04 1.0000E-03 3.0000E-02 0.0000E+00 p(1).n(1).p3
13 2.8815E+02 2.9312E+02 3.5815E+02 3.9285E+01 p(2).n(3).h2_temp
14 2.9315E+02 2.9315E+02 3.5815E+02 -7.2493E-01 p(2).n(3).wall_temp
13 2.8815E+02 2.9312E+02 3.5815E+02 3.9285E+01 p(2).n(3).h2_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
29 0.0000E+00 -7.5598E-04 0.0000E+00 7.5598E-04 p(2).n(4).Eqn(1): 0 = $h2_temp-(((((p2)*((1+alpha))))*((((temp_star-h2_temp))/((t_star+((time)*(p2))))))))
Variable Lower Value Upper $Value Name
0 -1.2346E+20 1.0000E+00 1.2346E+20 0.0000E+00 p(2).n(4).time
0 6.0000E+01 6.0000E+01 3.0000E+02 0.0000E+00 p(1).n(1).p2
0 5.0000E-04 1.0000E-03 3.0000E-02 0.0000E+00 p(1).n(1).p3
15 2.8815E+02 2.9320E+02 3.5815E+02 1.8882E+01 p(2).n(4).h2_temp
16 2.9315E+02 2.9315E+02 3.5815E+02 9.9172E-01 p(2).n(4).wall_temp
15 2.8815E+02 2.9320E+02 3.5815E+02 1.8882E+01 p(2).n(4).h2_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
2006 0.0000E+00 -1.0946E-01 0.0000E+00 1.0946E-01 p(1).c(2).t(2): not available
Variable Lower Value Upper $Value Name
2 2.9315E+02 2.9315E+02 3.5815E+02 -7.9554E+01 p(1).n(2).wall_temp
4 2.9315E+02 2.9315E+02 3.5815E+02 -4.1620E+01 p(1).n(3).wall_temp
6 2.9315E+02 2.9315E+02 3.5815E+02 -1.8591E+01 p(1).n(4).wall_temp
8 2.9315E+02 2.9315E+02 3.5815E+02 -8.8200E+00 p(1).n(5).wall_temp
10 2.9315E+02 2.9316E+02 3.5815E+02 -6.0293E+00 p(1).n(6).wall_temp
2 2.9315E+02 2.9315E+02 3.5815E+02 -7.9554E+01 p(1).n(2).wall_temp
0 2.9315E+02 2.9315E+02 3.5815E+02 0.0000E+00 p(1).n(1).wall_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
2007 0.0000E+00 -2.5022E-01 0.0000E+00 2.5022E-01 p(1).c(2).t(3): not available
Variable Lower Value Upper $Value Name
2 2.9315E+02 2.9315E+02 3.5815E+02 -7.9554E+01 p(1).n(2).wall_temp
4 2.9315E+02 2.9315E+02 3.5815E+02 -4.1620E+01 p(1).n(3).wall_temp
6 2.9315E+02 2.9315E+02 3.5815E+02 -1.8591E+01 p(1).n(4).wall_temp
8 2.9315E+02 2.9315E+02 3.5815E+02 -8.8200E+00 p(1).n(5).wall_temp
10 2.9315E+02 2.9316E+02 3.5815E+02 -6.0293E+00 p(1).n(6).wall_temp
4 2.9315E+02 2.9315E+02 3.5815E+02 -4.1620E+01 p(1).n(3).wall_temp
0 2.9315E+02 2.9315E+02 3.5815E+02 0.0000E+00 p(1).n(1).wall_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
2008 0.0000E+00 -3.3207E-01 0.0000E+00 3.3207E-01 p(1).c(2).t(4): not available
Variable Lower Value Upper $Value Name
2 2.9315E+02 2.9315E+02 3.5815E+02 -7.9554E+01 p(1).n(2).wall_temp
4 2.9315E+02 2.9315E+02 3.5815E+02 -4.1620E+01 p(1).n(3).wall_temp
6 2.9315E+02 2.9315E+02 3.5815E+02 -1.8591E+01 p(1).n(4).wall_temp
8 2.9315E+02 2.9315E+02 3.5815E+02 -8.8200E+00 p(1).n(5).wall_temp
10 2.9315E+02 2.9316E+02 3.5815E+02 -6.0293E+00 p(1).n(6).wall_temp
6 2.9315E+02 2.9315E+02 3.5815E+02 -1.8591E+01 p(1).n(4).wall_temp
0 2.9315E+02 2.9315E+02 3.5815E+02 0.0000E+00 p(1).n(1).wall_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
2009 0.0000E+00 -3.6373E-01 0.0000E+00 3.6373E-01 p(1).c(2).t(5): not available
Variable Lower Value Upper $Value Name
2 2.9315E+02 2.9315E+02 3.5815E+02 -7.9554E+01 p(1).n(2).wall_temp
4 2.9315E+02 2.9315E+02 3.5815E+02 -4.1620E+01 p(1).n(3).wall_temp
6 2.9315E+02 2.9315E+02 3.5815E+02 -1.8591E+01 p(1).n(4).wall_temp
8 2.9315E+02 2.9315E+02 3.5815E+02 -8.8200E+00 p(1).n(5).wall_temp
10 2.9315E+02 2.9316E+02 3.5815E+02 -6.0293E+00 p(1).n(6).wall_temp
8 2.9315E+02 2.9315E+02 3.5815E+02 -8.8200E+00 p(1).n(5).wall_temp
0 2.9315E+02 2.9315E+02 3.5815E+02 0.0000E+00 p(1).n(1).wall_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
2010 0.0000E+00 -3.8212E-01 0.0000E+00 3.8212E-01 p(1).c(2).t(6): not available
Variable Lower Value Upper $Value Name
2 2.9315E+02 2.9315E+02 3.5815E+02 -7.9554E+01 p(1).n(2).wall_temp
4 2.9315E+02 2.9315E+02 3.5815E+02 -4.1620E+01 p(1).n(3).wall_temp
6 2.9315E+02 2.9315E+02 3.5815E+02 -1.8591E+01 p(1).n(4).wall_temp
8 2.9315E+02 2.9315E+02 3.5815E+02 -8.8200E+00 p(1).n(5).wall_temp
10 2.9315E+02 2.9316E+02 3.5815E+02 -6.0293E+00 p(1).n(6).wall_temp
10 2.9315E+02 2.9316E+02 3.5815E+02 -6.0293E+00 p(1).n(6).wall_temp
0 2.9315E+02 2.9315E+02 3.5815E+02 0.0000E+00 p(1).n(1).wall_temp
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
2017 0.0000E+00 -7.1275E-04 0.0000E+00 7.1275E-04 p(2).c(2).t(3): not available
Variable Lower Value Upper $Value Name
12 2.9315E+02 2.9315E+02 3.5815E+02 -3.6833E+00 p(2).n(2).wall_temp
14 2.9315E+02 2.9315E+02 3.5815E+02 -7.2493E-01 p(2).n(3).wall_temp
16 2.9315E+02 2.9315E+02 3.5815E+02 9.9172E-01 p(2).n(4).wall_temp
18 2.9315E+02 2.9315E+02 3.5815E+02 1.6707E+00 p(2).n(5).wall_temp
20 2.9315E+02 2.9316E+02 3.5815E+02 1.8694E+00 p(2).n(6).wall_temp
14 2.9315E+02 2.9315E+02 3.5815E+02 -7.2493E-01 p(2).n(3).wall_temp
10 2.9315E+02 2.9316E+02 3.5815E+02 -6.0293E+00 p(1).n(6).wall_temp
************************************************
Also, this took around 2 minutes to run, so if there is any advice on how to decrease the computation time, that would be much appreciated!
You can see the model that Gekko writes with m.open_folder() and then open the model gk_model0.apm with a text editor:
Model
Parameters
p1 = 60.0, <= 300.0, >= 60.0
p2 = 0.001, <= 0.03, >= 0.0005
p3 = 0.17900000000000002, <= 1.611, >= 0.08950000000000001
p4 = 293.15, <= 303.15, >= 288.15
p5
End Parameters
Variables
v1 = 288.15, <= 358.15, >= 288.15
v2 = 293.15, <= 358.15, >= 293.15
v3 = 0.17900000000000002, <= 1.79, >= p3
v4 = 1000000.0, <= 35000000.0, >= 0.0
End Variables
Intermediates
i0=(1.0272572651140832/p2)
i1=(1.416731291198139/(1+i0))
i2=((i0)/((1+i0)))
i3=((p3)/(p2))
i4=2.762640041099948
i5=(((i1)*(p4))+((i2)*(v2)))
i6=(1+((0.02393942687999997)*((1-((((v1)/(33.14999999999998)))^(0.5))))))
i7=(0.074/v3)
End Intermediates
Equations
v3=((p1)*((p3+((p2)*([0. 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13
0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27
0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4 0.41
0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51 0.52 0.53 0.54 0.55
0.56 0.57 0.58 0.59 0.6 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69
0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8 0.81 0.82 0.83
0.84 0.85 0.86 0.87 0.88 0.89 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97
0.98 0.99 1. ])))))
$v1=((((p1)*((1+i0))))*((((i5-v1))/((i3+[0. 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13
0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27
0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4 0.41
0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51 0.52 0.53 0.54 0.55
0.56 0.57 0.58 0.59 0.6 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69
0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8 0.81 0.82 0.83
0.84 0.85 0.86 0.87 0.88 0.89 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97
0.98 0.99 1. ])))))
$v2=((p1)*((((v1-v2))/(i4))))
((v4)*(1000000.0))=((p1)*((((((4.1158415841584155)*(v1)))/((i7-8.165418118811874e-06)))-((((5.036651227998414e-09)*(((i6)^(2)))))/((((i7)*((i7+8.165418118811874e-06)))+((8.165418118811874e-06)*((i7-8.165418118811874e-06)))))))))
((v1)*(p5))<=358.15
((v2)*(p5))<=358.15
minimize p1
End Equations
Connections
p(100).n(end).v3=1.79
p(100).n(end).v3=fixed
p(100).n(end).v4=35000000.0
p(100).n(end).v4=fixed
End Connections
End Model
The problem is with the first two equation that have m.time and tm as Numpy arrays instead of using a Gekko variable or parameter. If a Numpy array or Python list is is needed in the optimization problem then create a new m.Param() such as:
t = m.Param(tm)
m.Equation(mass==t_f*(m_0 + m_dot_in*t))
m.Equation(T.dt()==t_f*(1 + alpha)*((T_star-T)/(t_star+t)))
Because the final time is minimized, the time in those equations may need to be t*t_f.
Equations such as m.Equation(T_w*final<=85+273.15) should be reformulated as m.Equation((T_w-(85+273.15))*final<=0) so that when final=0 then it is 0<=0. In this case, your original equations are okay but it is good practice to put all terms on one side of the equation.
Even with these modifications there is still Exception: @error: Solution Not Found. The problem may be with the terminal constraint. One way to get a feasible solution is to "soften" the constraint by including it as an objective.
m.Minimize(final*(mass-m_f)**2)
m.Minimize(final*(p-35.0e6)**2)
There is still a message that the problem is not feasible. You may want to continue to simplify your problem by turning off degrees of freedom (STATUS=0) and remove inequality constraints just to see if there is a problem with one of the equations such as divide-by-zero.
# HYDROGEN TANK REFUELING MODEL
# OPTIMIZE MODEL WITH GEKKO OPTIMIZATION SUITE
from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
M_H2 = 2.02
# CONSTANTS
## TANK PARAMETERS (ASSUME TYPE III, ALUMINUM, 74 l, RATED FOR 35 MPa)
V = 0.074 # m^3
a_in = 167/19e-3 # W/m^2/K
c_w = 2730 # J/kg/K
rho_w = 900 # kg/m^3
m_w = rho_w*(np.pi*(((0.358+19e-3)/2)**2)*(0.893+19e-3) - np.pi*((0.358/2)**2)*0.893)
A_in = 2*np.pi*(0.358/2)*((0.358/2) + 0.893) # m^2
T_w0 = 293.15 # K
m_f = 1.79 # final mass of hydrogen in tank, kg
## HYDROGEN PARAMETERS
c_p = 14.615e3 # specific heat at constant pressure, J/kg/K
c_v = 10.316e3 # specific heat at constant volume, J/kg/K
gamma = c_p/c_v
R = 8.314/M_H2 # gas constant for H2, J/kgK
T_c = -240 + 273.15 # critical temperature for H2, K
p_c = 1.3e6 # critical pressure for H2, Pa
w_H2 = -0.219 # acentric factor for H2
a = 0.45724*(R**2 * T_c**2)/(p_c**2)
b = 0.0778*(R*T_c)/p_c
kappa = 0.37464 + 1.54226*w_H2 - 0.26992*(w_H2**2)
# SET TIME ANALYSIS POINTS
nt = 101
tm = np.linspace(0, 1, nt)
m.time = tm
t = m.Param(tm)
# options
m.options.NODES = 6
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.MAX_ITER = 500
m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0
# FIXED VARIABLES
t_f = m.FV(value=60.0,lb=60.0,ub=300.0) # final fuel time, s
t_f.STATUS = 0
m_dot_in = m.FV(value=0.001,lb=0.0005,ub=0.03) # mass flow rate into tank, kg/s
m_dot_in.STATUS = 1
m_0 = m.FV(value=0.1*m_f,lb=0.05*m_f,ub=0.9*m_f) # initial mass of H2 in tank (as % of m_f), kg
m_0.STATUS = 1
T_inf = m.FV(value=20 + 273.15,lb=15 + 273.15,ub=30 + 273.15) # inflow temperature, K
T_inf.STATUS = 1
# PARAMETERS
f = np.zeros(nt)
f[-1] = 1.0
final = m.Param(value=f)
# VARIABLES
T = m.Var(value=15+273.15,lb=15+273.15,ub=85+273.15)
T_w = m.Var(value=T_w0,lb=T_w0,ub=85+273.15)
mass = m.Var(value=m_0,lb=m_0,ub=m_f)
p = m.Var(value=1.0e6,lb=0.0,ub=35.0e6)
# INTERMEDIATES
alpha = m.Intermediate((a_in*A_in)/c_v/m_dot_in)
gamma_p = m.Intermediate(gamma/(1 + alpha))
alpha_p = m.Intermediate(alpha/(1 + alpha))
t_star = m.Intermediate(m_0/m_dot_in)
t_w_star = m.Intermediate((m_w*c_w)/(a_in*A_in))
T_star = m.Intermediate(gamma_p*T_inf + alpha_p*T_w)
alpha_T = m.Intermediate(1 + kappa*(1 - (T/T_c)**0.5))
v = m.Intermediate(V/mass) # specific volume, m^3/kg
# EQUATIONS
m.Equation(mass==t_f*(m_0 + m_dot_in*t))
m.Equation(T.dt()==t_f*(1 + alpha)*((T_star-T)/(t_star+t)))
m.Equation(T_w.dt()==t_f*((T-T_w)/t_w_star))
m.Equation(p*1.0e6==t_f*((R*T/(v-b)) - ((a*alpha_T**2)/(v*(v+b) + b*(v-b)))))
m.Equation((T-(85+273.15))*final<=0)
m.Equation((T_w-(85+273.15))*final<=0)
# SPECIFIY ENDPOINT CONDITIONS
m.Minimize(final*(mass-m_f)**2)
m.Minimize(final*(p-35.0e6)**2)
#m.fix(mass, pos=len(m.time)-1, val=m_f)
#m.fix(p, pos=len(m.time)-1, val=35.0e6)
# MINIMIZE FINAL FUEL TIME
m.Minimize(t_f)
# SOLVE
m.options.SOLVER = 3
m.solve()
# RESULTS
print('Final Time: ' + str(t_f.value[0]))
Response to edit
The infeasibility reveals that wall_temp is a culprit from the first block. The solver is trying to push that value lower but it is at a bound. There are other parameters (variable 0) in the equation but wall_temp is the only one that is a variable that is at the lower bound. I made a modification to create a deadband to penalize any deviation outside of the upper (SPHI) and lower (SPLO) bounds. This way, the solution can still remain feasible if the constraints cannot be satisfied. The weights (WSPHI / WSPLO) can be increased if needed. Here is additional information on tuning.
# HYDROGEN TANK REFUELING MODEL
# OPTIMIZE MODEL WITH GEKKO OPTIMIZATION SUITE
from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
# CONSTANTS
## TANK PARAMETERS (ASSUME TYPE III, ALUMINUM, 74 l, RATED FOR 35 MPa)
V = 0.074 # m^3
a_in = 167/19e-3 # W/m^2/K
c_w = 2730 # J/kg/K
rho_w = 900 # kg/m^3
m_w = rho_w*(np.pi*(((0.358+19e-3)/2)**2)*(0.893+19e-3) - np.pi*((0.358/2)**2)*0.893)
A_in = 2*np.pi*(0.358/2)*((0.358/2) + 0.893) # m^2
T_w0 = 293.15 # K
m_f = 1.79 # final mass of hydrogen in tank, kg
## HYDROGEN PARAMETERS
c_p = 14.615e3 # specific heat at constant pressure, J/kg/K
c_v = 10.316e3 # specific heat at constant volume, J/kg/K
gamma = c_p/c_v
M_H2 = 2.02
R = 8.314/M_H2 # gas constant for H2, J/kgK
T_c = -240 + 273.15 # critical temperature for H2, K
p_c = 1.3e6 # critical pressure for H2, Pa
w_H2 = -0.219 # acentric factor for H2
a = 0.45724*(R**2 * T_c**2)/(p_c**2)
b = 0.0778*(R*T_c)/p_c
kappa = 0.37464 + 1.54226*w_H2 - 0.26992*(w_H2**2)
## SET INFLOW TEMPERATURE AND INITIAL MASS IN TANK
m_0 = 0.1*m_f
T_inf = 20 + 273.15
# SET TIME ANALYSIS POINTS
nt = 101
tm = np.linspace(0, 1, nt)
m.time = tm
t = m.Param(tm, name='time')
# options
m.options.NODES = 6
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.MAX_ITER = 500
m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0
# FIXED VARIABLES
t_f = m.FV(value=60.0,lb=60.0,ub=300.0) # final fuel time, s
t_f.STATUS = 1
m_dot_in = m.FV(value=0.001,lb=0.0005,ub=0.03) # mass flow rate into tank, kg/s
m_dot_in.STATUS = 1
# m_0 = m.FV(value=0.1*m_f,lb=0.05*m_f,ub=0.9*m_f) # initial mass of H2 in tank (as % of m_f), kg
# m_0.STATUS = 0
# T_inf = m.FV(value=20 + 273.15,lb=15 + 273.15,ub=30 + 273.15) # inflow temperature, K
# T_inf.STATUS = 0
# PARAMETERS
f = np.zeros(nt)
f[-1] = 1.0
final = m.Param(value=f, name='final')
# VARIABLES
T = m.CV(value=15+273.15,name='H2 Temp')
T.SPLO=15+273.15
T.SPHI=85+273.15
T.WSPLO = 100
T.WSPHI = 100
T.TR_INIT = 0
T.STATUS = 1
T_w = m.CV(value=T_w0,name='Wall Temp')
T_w.SPLO=T_w0
T_w.SPHI=85+273.15
T_w.WSPLO = 100
T_w.WSPHI = 100
T_w.TR_INIT = 0
T_w.STATUS = 1
# INTERMEDIATES
alpha = m.Intermediate((a_in*A_in)/c_v/m_dot_in, name='alpha')
gamma_p = m.Intermediate(gamma/(1 + alpha), name='gamma_p')
alpha_p = m.Intermediate(alpha/(1 + alpha), name='alpha_p')
t_star = m.Intermediate(m_0/m_dot_in, name='t_star')
t_w_star = m.Intermediate((m_w*c_w)/(a_in*A_in), name='t_w_star')
T_star = m.Intermediate(gamma_p*T_inf + alpha_p*T_w, name='Temp_star')
# alpha_T = m.Intermediate(1 + kappa*(1 - (T/T_c)**0.5))
# v = m.Intermediate(V/mass) # specific volume, m^3/kg
# EQUATIONS
# m.Equation(mass==t_f*(m_0 + m_dot_in*t*t_f))
m.Equation(T.dt()==t_f*(1 + alpha)*((T_star-T)/(t_star+t*t_f)))
m.Equation(T_w.dt()==t_f*((T-T_w)/t_w_star))
# m.Equation(p==t_f*((R*T/(v-b)) - ((a*alpha_T**2)/(v*(v+b) + b*(v-b)))))
# m.Equation((T-(85+273.15))*final<=0)
# m.Equation((T_w-(85+273.15))*final<=0)
# SPECIFIY ENDPOINT CONDITIONS
# m.Minimize(final*(mass-m_f)**2)
# m.Minimize(final*(p-35.0e6)**2)
m.Minimize(final*(T-351)**2)
#m.fix(mass, pos=len(m.time)-1, val=m_f)
#m.fix(p, pos=len(m.time)-1, val=35.0e6)
# MINIMIZE FINAL FUEL TIME
m.Minimize(t_f)
# SOLVE
m.options.SOLVER = 3
m.solve()
# RESULTS
print('Final Time: ' + str(t_f.value[0]))
Nice work simplifying the problem. Those deleted equations can likely be added back in, if desired, for convenience in plotting the solution. The next step is likely to plot the solution and see if it make intuitive sense with the weights and any constraint violations.