氢能汽车加注模型优化
Hydrogen Vehicle Refueling Model Optimization
我目前正在尝试确定使用 GEKKO 的氢 (H2) 汽车加氢过程的最佳流入条件。下面是耦合的常微分方程，它们控制 H2 和燃料箱壁的温度如何随加燃料时间而变化。
T.dt() = (1+alpha)*(T_star - T)/(t_star + t)
T_w.dt() = (T - T_w)/t_w_star
哪里
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)
这里，m_0是罐中H2的初始质量，m_dot_in是进入罐中的H2的质量流量，gamma是比热比对于H2，T_inf为H2的流入温度，其他变量为中间variables/tank参数。通过加油过程，m_dot_in 被取为常数（但未知），因此罐中 H2 的质量随时间的变化定义为：
m = m_0 + m_dot_in*t
此外，罐内 H2 的压力可以用真实的气体状态方程计算（我对这个模型使用 Peng-Robinson 状态方程）。
我试图用这个模型做的是确定最佳 m_dot_in、T_inf 和 m_0 以最小化总加油时间 t_f。对变量的一些限制是 T<=358.15 K 在整个加油过程中（出于安全原因），并且罐内 H2 的最终压力必须为 35 MPa。对于这个模型，我认为 t_f、m_dot_in、m_0 和 T_inf 是具有以下界限的固定变量：
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
下面我使用 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]))
此代码目前出现以下错误：
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...
总的来说，我对优化还很陌生，我尝试包括几个不同的等式和不等式约束，但似乎没有任何效果。我认为我根据 APMonitor 网站上的示例问题和信息正确地做了这件事，但显然我的实施有些问题。我想知道是否有人知道我应该做什么 change/add 或者我是否做错了什么？任何帮助将不胜感激！
感谢您的宝贵时间，
埃文
编辑：根据何登仁博士的回答，我尝试简化模型，不包括mass和p的变量，因为它们仅取决于 t_f、m_dot_in 和 T 的最终值，并且可以在获得解后计算。下面是我编辑的代码：
# 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]))
我仍然遇到不可行性（不像以前那么多），但我无法理解所说的不可行性的含义以及如何解决它们。以下是我得到的不可行性：
************************************************
***** 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
************************************************
此外，运行 大约需要 2 分钟，所以如果有任何关于如何减少计算时间的建议，我们将不胜感激！
可以看到Gekko用m.open_folder()写的模型，然后用文本编辑器打开模型gk_model0.apm：
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
问题出在前两个方程 m.time 和 tm 作为 Numpy 数组而不是使用 Gekko 变量或参数。如果在优化问题中需要 Numpy 数组或 Python 列表，则创建一个新的 m.Param() 例如：
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)))
因为最后的时间是最小化的，所以这些等式中的时间可能需要t*t_f。
m.Equation(T_w*final<=85+273.15) 等式应重新表述为 m.Equation((T_w-(85+273.15))*final<=0)，以便当 final=0 时 0<=0。在这种情况下，您的原始方程式没问题，但最好将所有项放在方程式的一侧。
即使进行了这些修改，仍然有 Exception: @error: Solution Not Found。问题可能出在终端约束上。获得可行解决方案的一种方法是通过将约束包含为 objective.
来“软化”约束
m.Minimize(final*(mass-m_f)**2)
m.Minimize(final*(p-35.0e6)**2)
还有提示问题不可行。您可能希望通过关闭自由度 (STATUS=0) 和删除不等式约束来继续简化您的问题，只是为了查看其中一个方程是否存在问题，例如被零除。
# 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]))
回复编辑
不可行性表明 wall_temp 是第一个区块的罪魁祸首。求解器正试图将该值推低，但它处于一定范围内。等式中还有其他参数（变量 0），但 wall_temp 是唯一一个处于下限的变量。我进行了修改以创建一个死区来惩罚上限 (SPHI) 和下限 (SPLO) 之外的任何偏差。这样，如果无法满足约束，解决方案仍然可行。如果需要，可以增加权重 (WSPHI / WSPLO)。这里是 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]))
很好地简化了问题。如果需要，可以重新添加那些删除的方程式，以便于绘制解决方案。下一步可能会绘制解决方案并查看它是否对权重和任何约束违规具有直观意义。
我目前正在尝试确定使用 GEKKO 的氢 (H2) 汽车加氢过程的最佳流入条件。下面是耦合的常微分方程，它们控制 H2 和燃料箱壁的温度如何随加燃料时间而变化。
T.dt() = (1+alpha)*(T_star - T)/(t_star + t)
T_w.dt() = (T - T_w)/t_w_star
哪里
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)
这里，m_0是罐中H2的初始质量，m_dot_in是进入罐中的H2的质量流量，gamma是比热比对于H2，T_inf为H2的流入温度，其他变量为中间variables/tank参数。通过加油过程，m_dot_in 被取为常数（但未知），因此罐中 H2 的质量随时间的变化定义为：
m = m_0 + m_dot_in*t
此外，罐内 H2 的压力可以用真实的气体状态方程计算（我对这个模型使用 Peng-Robinson 状态方程）。
我试图用这个模型做的是确定最佳 m_dot_in、T_inf 和 m_0 以最小化总加油时间 t_f。对变量的一些限制是 T<=358.15 K 在整个加油过程中（出于安全原因），并且罐内 H2 的最终压力必须为 35 MPa。对于这个模型，我认为 t_f、m_dot_in、m_0 和 T_inf 是具有以下界限的固定变量：
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
下面我使用 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]))
此代码目前出现以下错误：
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...
总的来说，我对优化还很陌生，我尝试包括几个不同的等式和不等式约束，但似乎没有任何效果。我认为我根据 APMonitor 网站上的示例问题和信息正确地做了这件事，但显然我的实施有些问题。我想知道是否有人知道我应该做什么 change/add 或者我是否做错了什么？任何帮助将不胜感激！
感谢您的宝贵时间，
埃文
编辑：根据何登仁博士的回答，我尝试简化模型，不包括mass和p的变量，因为它们仅取决于 t_f、m_dot_in 和 T 的最终值，并且可以在获得解后计算。下面是我编辑的代码：
# 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]))
我仍然遇到不可行性（不像以前那么多），但我无法理解所说的不可行性的含义以及如何解决它们。以下是我得到的不可行性：
************************************************
***** 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
************************************************
此外，运行 大约需要 2 分钟，所以如果有任何关于如何减少计算时间的建议，我们将不胜感激！
可以看到Gekko用m.open_folder()写的模型，然后用文本编辑器打开模型gk_model0.apm：
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
问题出在前两个方程 m.time 和 tm 作为 Numpy 数组而不是使用 Gekko 变量或参数。如果在优化问题中需要 Numpy 数组或 Python 列表，则创建一个新的 m.Param() 例如：
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)))
因为最后的时间是最小化的，所以这些等式中的时间可能需要t*t_f。
m.Equation(T_w*final<=85+273.15) 等式应重新表述为 m.Equation((T_w-(85+273.15))*final<=0)，以便当 final=0 时 0<=0。在这种情况下，您的原始方程式没问题，但最好将所有项放在方程式的一侧。
即使进行了这些修改，仍然有 Exception: @error: Solution Not Found。问题可能出在终端约束上。获得可行解决方案的一种方法是通过将约束包含为 objective.
m.Minimize(final*(mass-m_f)**2)
m.Minimize(final*(p-35.0e6)**2)
还有提示问题不可行。您可能希望通过关闭自由度 (STATUS=0) 和删除不等式约束来继续简化您的问题，只是为了查看其中一个方程是否存在问题，例如被零除。
# 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]))
回复编辑
不可行性表明 wall_temp 是第一个区块的罪魁祸首。求解器正试图将该值推低，但它处于一定范围内。等式中还有其他参数（变量 0），但 wall_temp 是唯一一个处于下限的变量。我进行了修改以创建一个死区来惩罚上限 (SPHI) 和下限 (SPLO) 之外的任何偏差。这样，如果无法满足约束，解决方案仍然可行。如果需要，可以增加权重 (WSPHI / WSPLO)。这里是 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]))
很好地简化了问题。如果需要，可以重新添加那些删除的方程式，以便于绘制解决方案。下一步可能会绘制解决方案并查看它是否对权重和任何约束违规具有直观意义。