问题：
氢汽车加油模型优化
缪风史
我目前正在尝试确定使用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的流入温度，其他变量是中间变量/储罐参数。通过加油过程，m_dot_in被认为是恒定的（但未知），因此储罐中H2随时间的质量定义为：
m = m_0 + m_dot_in*t
此外，罐内H2的压力可以用真实的气体状态方程来计算（我在这个模型中使用彭-罗宾逊状态方程）。
我试图用这个模型做的是确定最佳的m_dot_in、T_inf和m_0以最小化总加油时间，t_f。变量的一些约束是
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网站上的示例问题和信息正确地完成的，但很明显我的实现出了问题。我想知道是否有人知道我应该更改/添加什么，或者我是否做了一些完全错误的事情？任何帮助都将不胜感激！
谢谢你抽出时间，
埃文
编辑：根据何登仁博士的回答，我试图简化模型，使质量和p的变量不包括在内，因为它们只取决于t\u f、m\u dot\u 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分钟的运行时间，因此如果有任何关于如何减少计算时间的建议，将不胜感激！
共有1个答案
尚俊楠
您可以看到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方程（T\u w*final
即使进行了这些修改，仍然存在异常：@错误：未找到解决方案。问题可能与终端约束有关。获得可行解的一种方法是通过将约束作为目标来“软化”约束。
m.Minimize(final*(mass-m_f)**2)
m.Minimize(final*(p-35.0e6)**2)
仍然有一个信息表明这个问题不可行。您可能希望通过关闭自由度（状态=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）。这是关于调整的其他信息。
# 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]))
简化问题的好工作。如果需要，可以将这些删除的方程添加回，以方便绘制解决方案。下一步可能会绘制解决方案，看看它是否对权重和任何违反约束的行为有直观的意义。
类似资料：
-
我实际上是在用R建立一个电动汽车时间偏离的预测模型，作为数据，我有两年前的时间偏离历史，我希望建立一个以日期（日）为输入的预测模型，模型的输出将是时间偏离。一个步骤可以考虑在两个开始之间间隔一个小时（例如:-->) > 我试图将问题看作是以（1表示“有一个离开”，0表示“无”）聚类一个，并应用随机森林模型（日和小时作为输入，1或0作为输出），但该模型找不到输出和输入之间的联系，而不考虑日或小时，该
-
优享型用车 基本描述 优享是滴滴出行为客户提供的更加舒适、经济的标准化出行服务。优享车型提供大空间车型，配备全认证高星级司机，旨在为企业提供更宽敞舒适的用车，更有品质的乘车体验。 车型优势 相比较于普通快车、专车，优享用车具有以下优势： 比快车更舒适：更严格的车价及轴距准入要求； 比专车更经济：亲民定价类似快车的计价结构； 标准化服务：专业的司机培训带来优质服务； 业务流程示例 接入流程 优享型用
-
电话面，周一上午给我打电话约的周二下午 面试官拷打的挺多的，我还是第一次遇到这么不温柔的面试官 问题一个接着一个的抛 我太菜了给我问蒙了都 对深蓝汽车了解吗 接受工作地点在重庆吗 ==和equals的区别 Integer a=200和Integer b=200用equals比较和用==比较分别是true还是false等等一系列（接连拷打） 迭代器模式，为啥不用for循环 工厂模式和抽象工厂模式 r
-
10-8 面试官懂的太多了，秋招以来第一次面试被问麻了。 手撕 简单dp题，到右下角的最短路径 面试官说可以把边界条件拿出来做，这样会更清晰点 Pytorch DDP了解过吗 不了解 CV的发展路径 从AlexNet开始说，因为想不起来具体改进，就总结了说是各种架构和激活函数的改进 NLP的发展路径 RNN-》LSTM-》Transformer 不清楚是不是这个发展 RNN和Transformer
-
上正文： 1. == equals (√) 2.出题：Integer a = 100,Integer b = 100.....balabala.....就是一些Integer和int这些的比较，equals，==(答上来了，但是他问我为啥。我说是equals比较的是值，他问我为啥。说实话我很懵逼。下来复盘的时候才get应该是想问我拆箱装箱，Integercache...其实是想到了但是没敢说) 3
-
#长安汽车2023秋招面经#一面-综合面-10-15min 自我介绍-优势、性格、爱好 为什么选择长安，对长安的理解 除了专业书还看什么书 周末怎么安排 身边的同学怎么评价你 在压力环境怎么处理工作 二面-10-15min 线下专业面试 因为是信息安全岗位，三个简历项目和安全相关所以就着重介绍一个和汽车会相关的安全项目，面试官问了一些防御的措施-不会但基本往会的说哈哈 之后问了有没有男朋友啥的（笑