问题描述
我目前正在使用 Gekko 库研究动态传热模型。该模型由偏微分方程组组成。如APMonitor平台所示，空间变化的微分已被离散化，而时间上的微分则通过有限元正交搭配法进行数值求解。
我有两个问题：
-
该模型仅使用空间中最多 3 个离散点和
NODES = 3时间中最多 30 个离散点求解。而对于NODES = 4，节点不遵循函数的值。什么解释了这种行为？ -
当使用任何 if2/3 时，模型都无法求解。测试了所有类型的求解器，但没有找到解决方案。当 m.if3 被 m.Intermediates 替换时，模型正常求解。此外，如果未指定节点数，则最多使用 n = 20 和 nt = 300 个点来求解模型。通过添加 if 和节点可以解释这种性能下降的原因是什么？模型的表述或结构是否不正确？
可以做些什么来改进？该模型需要至少用 100 个时间点和 NODES = 3 求解。
***我也想知道如何在另一个 if2/3 中使用 if2/3。
这是模型的主要部分：
import numpy as np
import data
from func_H import func_h
from func_KW import func_kw
from func_den import func_fg_cpv
from func_PG import func_pg
from func_CPV import func_cpv11,func_cpv21,func_cpv31,func_cpv41,func_cpv51,func_cpv6,func_cpv7,func_cpv8,func_cpv12,func_cpv22,func_cpv32,func_cpv42,func_cpv52
from func_CPL import func_cpl11,func_cpl12,func_cpl21,func_cpl22,func_cpl31,func_cpl32,func_cpl4,func_cpl5,func_cpl6,func_cpl71,func_cpl72,func_cpl8,func_cpl9
from func_cond_in import func_cond_inic
from func_DHW import func_dhw
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=False)
zf = 0.17
tf = 768.24
dz = zf/n
n = 3
nt = 10
m.time = np.linspace(0,tf,nt)
nd = 3
####################################################################################################
u = func_cond_inic(data.tg,data.g,data.p1,data.yi[0],data.yi[1],data.yi[2],data.yi[3],data.yi[4],data.yi[5],data.yi[6],data.yi[7],n)
tg = [m.Var(u[0][i]) for i in range(n)]
g = [m.Var(u[1][i]) for i in range(n)]
pg = [m.Var(u[2][i]) for i in range(n)]
tp = [m.Var(data.tp) for i in range(n)]
y1 = [m.Var(u[3][i][0]) for i in range(n)]
y2 = [m.Var(u[4][i][0]) for i in range(n)]
y3 = [m.Var(u[5][i][0]) for i in range(n)]
y4 = [m.Var(u[6][i][0]) for i in range(n)]
y5 = [m.Var(u[7][i][0]) for i in range(n)]
y6 = [m.Var(u[8][i][0]) for i in range(n)]
y7 = [m.Var(u[9][i][0]) for i in range(n)]
y8 = [m.Var(u[10][i][0]) for i in range(n)]
y = [np.array([y1[i],y2[i],y3[i],y4[i],y5[i],y6[i],y7[i],y8[i]]) for i in range(n)]
yt = [np.transpose(y[i]) for i in range(n)]
x1 = [m.Var(data.xi[0]) for i in range(n)]
x2 = [m.Var(data.xi[1]) for i in range(n)]
x3 = [m.Var(data.xi[2]) for i in range(n)]
x4 = [m.Var(data.xi[3]) for i in range(n)]
x5 = [m.Var(data.xi[4]) for i in range(n)]
x6 = [m.Var(data.xi[5]) for i in range(n)]
x7 = [m.Var(data.xi[6]) for i in range(n)]
x8 = [m.Var(data.xi[7]) for i in range(n)]
x9 = [m.Var(data.xi[8]) for i in range(n)]
x = [np.array([x1[i],x2[i],x3[i],x4[i],x5[i],x6[i],x7[i],x8[i],x9[i]]) for i in range(n)]
####################################################################################################
rhog = [m.Intermediate((pg[i]*(data.mmg @ y[i])/(data.r*tg[i]))) for i in range(n)]
kw = [m.Intermediate(func_kw(g[i],tg[i],y[i],pg[i],rhog[i])) for i in range(n)]
weg = [m.Intermediate(m.exp(-3.5+(0.0549/tg[i]))) for i in range(n)]
wg = [m.Intermediate((pg[i]*data.mmg[3]*y4[i])/(data.r*tg[i])) for i in range(n)]
rw_f = [m.Intermediate(kw[i]*data.a*(weg[i]-wg[i])) for i in range(n)]
rw_m = [m.Intermediate(rw_f[i]*1000/18.01528) for i in range(n)]
dhw = [m.Intermediate(func_dhw(tg[i],tp[i])) for i in range(n)]
qcond = [m.Intermediate(rw_m[i]*dhw[i]) for i in range(n)]
# # Heat capacity
cpv1 = [m.if2(tg[i]-700,func_cpv11(tg[i]),func_cpv12(tg[i])) for i in range(n)]
cpv2 = [m.if2(tg[i]-500,func_cpv21(tg[i]),func_cpv22(tg[i])) for i in range(n)]
cpv3 = [m.if2(tg[i]-1200,func_cpv31(tg[i]),func_cpv32(tg[i])) for i in range(n)]
cpv4 = [m.if2(tg[i]-1700,func_cpv41(tg[i]),func_cpv42(tg[i])) for i in range(n)]
cpv5 = [m.if2(tg[i]-1300,func_cpv51(tg[i]),func_cpv52(tg[i])) for i in range(n)]
# cpv1 = [m.Intermediate(func_cpv11(tg[i])) for i in range(n)]
# cpv2 = [m.Intermediate(func_cpv21(tg[i])) for i in range(n)]
# cpv3 = [m.Intermediate(func_cpv31(tg[i])) for i in range(n)]
# cpv4 = [m.Intermediate(func_cpv41(tg[i])) for i in range(n)]
# cpv5 = [m.Intermediate(func_cpv51(tg[i])) for i in range(n)]
cpv6 = [m.Intermediate(func_cpv6(tg[i])) for i in range(n)]
cpv7 = [m.Intermediate(func_cpv7(tg[i])) for i in range(n)]
cpv8 = [m.Intermediate(func_cpv8(tg[i])) for i in range(n)]
cpv = [[cpv1[i],cpv2[i],cpv3[i],cpv4[i],cpv5[i],cpv6[i],cpv7[i],cpv8[i]] for i in range(n)]
fg_cpv = [m.Intermediate(func_fg_cpv(g[i],cpv[i],y[i])) for i in range(n)]
h = [m.Intermediate(func_h(g[i],cpv[i])) for i in range(n)]
qconv = [m.Intermediate((data.a*h[i]*(tg[i]-tp[i]))) for i in range(n)]
p = [m.Intermediate(func_pg(g[i],rhog[i])) for i in range(n)]
# # Heat capacity
cpl1 = [m.if2(tp[i]-950,func_cpl11(tp[i]),func_cpl12(tp[i])) for i in range(n)]
cpl2 = [m.if2(tp[i]-900,func_cpl21(tp[i]),func_cpl22(tp[i])) for i in range(n)]
cpl3 = [m.if2(tp[i]-847,func_cpl31(tp[i]),func_cpl32(tp[i])) for i in range(n)]
# cpl1 = [m.Intermediate(func_cpl11(tp[i])) for i in range(n)]
# cpl2 = [m.Intermediate(func_cpl21(tp[i])) for i in range(n)]
# cpl3 = [m.Intermediate(func_cpl31(tp[i])) for i in range(n)]
cpl4 = [m.Intermediate(func_cpl4(tp[i])) for i in range(n)]
cpl5 = [m.Intermediate(func_cpl5(tp[i])) for i in range(n)]
cpl6 = [m.Intermediate(func_cpl6(tp[i])) for i in range(n)]
cpl7 = [m.if2(tp[i]-500,func_cpl71(tp[i]),func_cpl72(tp[i])) for i in range(n)]
# cpl7 = [m.Intermediate(func_cpl71(tp[i])) for i in range(n)]
cpl8 = [m.Intermediate(func_cpl8(tp[i])) for i in range(n)]
cpl9 = [m.Intermediate(func_cpl9(tp[i])) for i in range(n)]
cpl = [[cpl1[i],cpl2[i],cpl3[i],cpl4[i],cpl5[i],cpl6[i],cpl7[i],cpl8[i],cpl9[i]] for i in range(n)]
rhop1 = [m.Intermediate((x1[i]*data.rhoc[0]*(1-data.e))) for i in range(n)]
rhop2 = [m.Intermediate((x2[i]*data.rhoc[1]*(1-data.e))) for i in range(n)]
rhop3 = [m.Intermediate((x3[i]*data.rhoc[2]*(1-data.e))) for i in range(n)]
rhop4 = [m.Intermediate((x4[i]*data.rhoc[3]*(1-data.e))) for i in range(n)]
rhop5 = [m.Intermediate((x5[i]*data.rhoc[4]*(1-data.e))) for i in range(n)]
rhop6 = [m.Intermediate((x6[i]*data.rhoc[5]*(1-data.e))) for i in range(n)]
rhop7 = [m.Intermediate((x7[i]*data.rhoc[6]*(1-data.e))) for i in range(n)]
rhop8 = [m.Intermediate((x8[i]*data.rhoc[7]*(1-data.e))) for i in range(n)]
rhop9 = [m.Intermediate((x9[i]*data.rhoc[8]*(1-data.e))) for i in range(n)]
rhop = [np.array([rhop1[i],rhop2[i],rhop3[i],rhop4[i],rhop5[i],rhop6[i],rhop7[i],rhop8[i],rhop9[i]]) for i in range(n)]
rho_cpl = [m.Intermediate((rhop[i] @ cpl[i])) for i in range(n)]
####################################################################################################
m.Equation(tg[0] == data.tg)
m.Equation([tg[i] == tg[i-1] + dz*(qcond[i-1]-qconv[i-1])/fg_cpv[i-1] for i in range(1,n)])
m.Equation(g[0] == data.g)
m.Equation([g[i] == g[i-1] + dz*rw_f[i-1] for i in range(1,n)])
m.Equation(pg[0] == data.p1)
m.Equation([pg[i] == pg[i-1] + dz*p[i-1] for i in range(1,n)])
m.Equation(y1[0] == data.yi[0])
m.Equation(y2[0] == data.yi[1])
m.Equation(y3[0] == data.yi[2])
m.Equation(y4[0] == data.yi[3])
m.Equation(y5[0] == data.yi[4])
m.Equation(y6[0] == data.yi[5])
m.Equation(y7[0] == data.yi[6])
m.Equation(y8[0] == data.yi[7])
m.Equation([y1[i] == y1[i-1] - dz*rw_f[i-1]*y1[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y2[i] == y2[i-1] - dz*rw_f[i-1]*y2[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y3[i] == y3[i-1] - dz*rw_f[i-1]*y3[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y4[i] == y4[i-1] + dz*rw_f[i-1]*(1-y4[i-1])/g[i-1] for i in range(1,n)])
m.Equation([y5[i] == y5[i-1] - dz*rw_f[i-1]*y5[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y6[i] == y6[i-1] - dz*rw_f[i-1]*y6[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y7[i] == y7[i-1] - dz*rw_f[i-1]*y7[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y8[i] == y8[i-1] - dz*rw_f[i-1]*y8[i-1]/g[i-1] for i in range(1,n)])
m.Equation([tp[i].dt() == ((qconv[i]-qcond[i])/rho_cpl[i]) for i in range(n)])
m.Equation([x1[i].dt() == 0 for i in range(n)])
m.Equation([x2[i].dt() == 0 for i in range(n)])
m.Equation([x3[i].dt() == 0 for i in range(n)])
m.Equation([x4[i].dt() == 0 for i in range(n)])
m.Equation([x5[i].dt() == 0 for i in range(n)])
m.Equation([x6[i].dt() == 0 for i in range(n)])
m.Equation([x7[i].dt() == (-rw_f[i]/(rhop1[i]+rhop2[i]+rhop3[i]+rhop4[i]+rhop5[i]+rhop6[i]+rhop7[i]+rhop8[i]+rhop9[i])) for i in range(n)])
m.Equation([x8[i].dt() == 0 for i in range(n)])
m.Equation([x9[i].dt() == 0 for i in range(n)])
####################################################################################################
m.options.imode = 4
m.options.soLVER = 1
m.options.NODES = nd
m.options.REDUCE = 2
m.options.MAX_ITER = 250
m.solve()
完整模型可从以下链接获得：
https://drive.google.com/drive/folders/1_AFnBPZ_mcEhwZb80a3VYOcg8PVXP5Qj?usp=sharing
提前致谢。
解决方法
对于复杂的问题，第一个策略是仅用小步长（例如 0.01）解决一个时间步长（2 个点）。使用 m.options.SOLVER=1 切换到 APOPT 求解器也有助于找到问题的根源，因为它能够更好地识别不可行性。
m.time = np.linspace(0,0.01,2)
求解器报告一个不成功的解，因为在一次迭代后目标是 Inf（无穷大）。
Number of state variables: 510
Number of total equations: - 402
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 108
----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------
Iter Objective Convergence
0 1.67459E+11 1.01245E+05
1 +Inf 1.01245E+05
No feasible solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 1.1939000000000002 sec
Objective : 108000.
Unsuccessful with error code 0
---------------------------------------------------
Creating file: infeasibilities.txt
Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
@error: Solution Not Found
运行目录显示附加信息。运行目录为 m.path 或在解决命令之前使用 m.open_folder() 打开。
m.open_folder()
m.solve()
打开 infeasibilities.txt 文件显示几个中间体被零除。
rhog = [m.Intermediate((pg[i]*(data.mmg @ y[i])/(data.r*tg[i]))) for i in range(n)]
wg = [m.Intermediate((pg[i]*data.mmg[3]*y4[i])/(data.r*tg[i])) for i in range(n)]
下限有助于避免被零除，例如：
tg = [m.Var(u[0][i],lb=0.1) for i in range(n)]
即使有这个下界，还有另一个表达式可以达到无穷大。这个策略可以用来连续发现模型的问题。