文档建议反馈控制台
首页
学习
活动
专区
圈层
工具
MCP广场
文章/答案/技术大牛
发布
社区首页 >问答首页 >氢汽车加油模型优化

问氢汽车加油模型优化
EN

Stack Overflow用户
提问于 2021-06-22 22:41:08
回答 1查看 104关注 0票数 0

我目前正在努力确定使用GEKKO的氢(H2)汽车加油过程的最佳入流条件。下面是耦合的常微分方程，控制H2和燃料箱壁的温度随加油时间的变化。

代码语言：javascript
复制
T.dt() = (1+alpha)*(T_star - T)/(t_star + t)
T_w.dt() = (T - T_w)/t_w_star

哪里

代码语言：javascript
复制
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质量被定义为：

代码语言：javascript
复制
m = m_0 + m_dot_in*t

另外，可以用真实的气体状态方程来计算储罐内H2的压力(本模型采用彭-鲁滨逊状态方程)。

我试图用这个模型来确定最优的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是具有以下边界的固定变量：

代码语言：javascript
复制
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复制了这个优化问题的代码：

代码语言：javascript
复制
# 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]))

此代码当前给出了以下错误：

代码语言：javascript
复制
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网站上的示例问题和信息，我做得很正确，但是很明显，我的实现是有问题的。我想知道是否有人知道我应该改变什么/添加什么，或者我做了什么完全错误的事情？任何帮助都将不胜感激！

谢谢您抽时间见我,

埃文

编辑：基于Hedengren博士的回答，我试图简化模型，使mass和p的变量不包括在内，因为它们只取决于t_f、m_dot_in和T的最终值，并且可以在得到解之后进行计算。下面是我编辑的代码：

代码语言：javascript
复制
# 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]))

我仍然不可能(不像以前那么多)，但是我很难理解所谓的不可行意味着什么，以及如何去修复它们。以下是我所得到的不可能：

代码语言：javascript
复制
************************************************
***** 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分钟的运行，所以如果有任何关于如何减少计算时间的建议，这将是非常感谢！

python
optimization
gekko
EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2021-06-23 18:03:37

您可以看到Gekko用m.open_folder()编写的模型，然后用文本编辑器打开模型gk_model0.apm：

代码语言：javascript
复制
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()，如：

代码语言：javascript
复制
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。问题可能是终端约束。获得可行解决方案的一种方法是将约束“软化”，将其作为一个目标。

代码语言：javascript
复制
m.Minimize(final*(mass-m_f)**2)
m.Minimize(final*(p-35.0e6)**2)

仍然有一个信息表明，这个问题是不可行的。你可能想继续简化你的问题，关闭自由度(STATUS=0)，消除不等式约束，只是看看是否有一个方程，如除以零的问题。

代码语言：javascript
复制
# 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)和下限(人运)界限。这样，如果不能满足这些限制，解决方案仍然是可行的。如果需要，可以增加权重(WSPHI / WSPLO)。这是关于调优的其他信息。

代码语言：javascript
复制
# 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]))

简化这个问题做得很好。如果愿意的话，这些被删除的方程可能会被重新加入，以便于绘制解决方案。下一步很可能绘制出解决方案，看看它是否对权重和任何违反约束的行为具有直观的意义。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：

https://stackoverflow.com/questions/68091462

复制
相关文章

相似问题

领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 粤公网安备44030502008569号

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档